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ABSTRACT 

This   note    describes    the   application  of    the    Carnahan-Starl ing-DeSantis   equation 
of    state    to  halogenated  hydrocarbon  refrigerants  and   their  mixtures.      A  com- 
plete and  consistent   set   of   thermodynamic  functions   is   derived  from    the    p-V-T 
equation  of    state    and    the   perfect    (ideal)    gas  heat   capacities.      A  thorough 
discussion  of  reference    states   is   included  for   both  pure  materials  and   their 
mixtures.      Although   this  model    exhibits  a   critical    point,    it   does  not  quanti- 
tatively  represent   properties   in  the    critical    region.      Despite    this  limita- 
tion,   this  model    can  represent  both  liquid  and  gaseous  mixtures  away   from 
their   own  critical   points,    even  at   conditions   near   to  and  above    the    critical 
points   of    their   components. 

Algorithms   and  FORTRAN  routines  for    the   use    of    this  model    are   presented  along 
with  the  numerical    coefficients  for  11   pure  refrigerants  and  7    mixtures. 
Routines   for    evaluating   the    coefficients   from    saturation   data    are    included. 
Several   examples  of   the   application  of    this   equation  of    state   are    presented   to 
demonstrate    its  versatility.      It    is    shown   to   predict      the      properties   of    pure 
materials  well   and  to  describe   the   detailed  features   of   mixtures,    both  phase 
diagrams   and   thermodynamic   properties.      The   average    deviation  from    the    tabu- 
lated  saturation   properties   of    the   11   pure    refrigerants    is   0.54  %    for    pres- 
sures,   0.09  %    for    liquid  volumes,     and   0.50  %    for   vapor   volumes. 
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"P'        P 


C  ,     (C  )  molar   heat   capacity    at   constant   pressure    (for   the   perfect   gas 


state) 
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Cq,    ci,    c«  coefficients  for    the   temperature   dependence    of   C     (defined  by 
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p*  arbitrary  reference    pressure 

p°  pressure   of   perfect    gas    state   at    specified  temperature  and 

volume 

l>low'    p  lower   and  upper  bound   on  pressure   for   a  physically  meaningful 

solution  of   equation  of    state 


R 
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molar   entropy 
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temperature 
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1  function  defined  by  Eqn.    6.1 

6,   A  change    in  a  quantity 

6 '  x'  ith   guess   for   iteration  variable 
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p..  chemical   potential    of   component    i    in  phase   j 

a  molecular   diameter 
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1.       INTRODUCTION 

This  Technical    Note   describes   the  application  of    an  equation  of    state   based  on 
the  hard-sphere   fluid  to  an   important    class   of    industrial    fluids,    the  halogen- 
ated  hydrocarbon  refrigerants  and  their  mixtures.      These    compounds,    including 
several    azeotropic  mixtures,   have   been   developed   over    the    past    several    decades 
and  are   now    the  dominant  working  fluids   in  refrigeration  and  heat  pumping 
equipment.      Recently,    considerable   interest  has   arisen  in   the   use    of   non- 
azeotropic  mixtures   as  working  fluids  because   of  potential    efficiency   improve- 
ments and  operational   advantages.      The  lack  of   precise   information  on  the 
thermophysical    behavior  of   these  mixtures  has  hindered  their  commercial    appli- 
cation.     There   is  a   need  for   a   general    property  formulation  that  would  be 
applicable   to  the  analysis  of   a  variety   of   pure  refrigerants  and  their 
mixtures. 

The  purpose    of    this  Note   is  to   show    how    an  equation  of    state  founded  on  a 
realistic  physical  model   can  be  used  to   describe    the    thermodynamic   properties 
of   a  pure  fluid  or  mixture  with  a  minimum   amount  of   experimental    information. 
Such  an  equation   can  be  viewed  on   two  levels.      First,    it  represents  a    correla- 
tion of    experimental    data.      But   more   than  merely   correlating   isolated  data,    an 
equation  of   state  models  the    thermophysical   behavior   of   a   fluid.      Such   a  model 
has  built   into   it    the  virtues  of    thermodynamic  consistency,    realistic  limiting 
behavior,    and  the  power  to  predict  unmeasured  properties  reasonably  when  there 
is  limited  experimental    information.      Unlike    simple  correlating  schemes,    an 
equation  of   state  founded  on  a  realistic   physical   model    becomes  an  increas- 
ingly powerful    tool    as   ones  knowledge    about   a  material    increases.      The  authors 


wish   to  emphasize,    however,    that   there   is  no   substitute   for   data*,  any 
predictive    scheme,    even  one   founded  on  a   good  physical   model,    can  be 
considered  only   a   temporary   substitute  for   experimental    information. 

This  Note   develops  a   particular   equation  of    state — a  modification  of    the 
Carnahan-Starl  ing  hard   sphere  fluid  that  was  first   proposed  by  DeSantis,    et 
al.    (1976): — and  applies  it  to  the  representation  of  refrigerant    thermodynamic 
properties.       We  will   refer    to   this  as   the  Carnahan-Starl  ing-DeSantis   (CSD) 
equation  of    state.      Much  of   this   development   is   general,    and  can  be  applied  to 
other    similar   equations  of    state.      While   only   refrigerants  are  discussed  here, 
the   CSD  equation  of   state  has  been  applied  to  other   classes  of   fluids,    such  as 
hydrocarbons   and    simple    inorganic   molecules    (DeSantis,    et   al.    1976). 

Section  2   of    this  Note   briefly   describes  and  compares   classes  of   equations  of 
state.      The  motivation  for    choosing  the  hard   sphere  reference    fluid  is 
explained.      Section  3   presents  mixing  rules  for    the   extension  of    the   equation 
of    state    to  mixtures  and  introduces  an  empirical    interaction  parameter  used  to 
account  for   deviations  from   the  Lorentz   mixing  rules.      Methods   of   predicting 
the   interaction  parameter  from  molecular  data  are  reviewed  and  the  need  for 
experimental   p-V-T-x  data  for  its   evaluation  is   emphasized.      A  complete    set   of 
thermodynamic  functions   is  developed   in  Section  4.      These   formulations  use 
only   the   p-V-T  information  of    the   equation  of    state  and   the  heat   capacity   of  a 
perfect    (ideal)    gas.      A  discussion  of    reference    states   is   also    included. 
Section  5   considers   the   behavior   of    the   CSD  and  other   equations   of    state   near 
the   critical    point.      In  Section  6,    the  numerical    implementation  of    the   equa- 
tion of   state  is  discussed,    including   the   determination  of    the   pure    component 


parameters  and  the  mixture   interaction  coefficient.      A  set  of   computer 
subroutines  for   carrying  out    these   calculations  is   presented.      Section  7 
examines   the  ability   of    the  equation  of   state   to  represent  the  thermodynamic 
characteristics   of  halogenated  hydrocarbon  refrigerants  and    their  mixtures. 
Concluding  remarks   are   presented  in  Section  8.      Finally,    the  Appendices   con- 
tain a   tabulation  of    coefficients  for   use  with   the  hard  sphere   equation  of 
state   and  a   listing  of    the  FORTRAN    subroutines   developed  for    its 
implementation. 


2.      A  COMPARISON  OF  EQUATIONS  OF  STATE 

An  equation  of    state  represents  the  pressure— volume-temperature  behavior  of   a 
pure  material   or  the  p-V-T-x  behavior  of  a  mixture.      When  combined  with  heat 
capacity    information  the  complete   thermodynamic  behavior   of  a  material 
(including   enthalpy,    entropy,    Gibbs   free   energy,    etc.)    is    defined.      For    the 
correlation  of    experimental    data,    an  equation  of    state    can  represent   all    the 
thermodynamic  properties  with   strict   consistency,    an  advantage   not   necessarily 
retained  when   separate   functions  are   used  to   correlate    individual    properties. 
An  equation  of    state   also   permits  reasonable   evaluation  of   properties  for 
which   there  are  no  measurements. 

There  are  a   number   of    classes  of    equations  of    state  which  are  based  on  varying 
combinations  of   theory  and  empiricism  and  with  varying  ranges   of 
applicability.      The   simplest   of    all   equations  of    state    is   the  perfect   gas  law: 


£V  =    i  (2.1) 

RT 


It   describes  the  behavior    of   an  assembly  of   point  masses  with  no  intermol  ecu- 
lar  attraction  or  repulsion.       It  represents   the   properties   of  a   gas   in   the 
limit   of   zero  pressure.      Depending  upon  ones  requirements,    it   is   useful    at 
pressures  less   than  0.1   MPa*,  at   temperatures  well   above    the    critical    tempera- 
ture   it    can  be    useful    at  higher   pressures. 

The  elegant  form  of   the   perfect    gas  law    is  retained  in  the  virial    equation  of 
state: 


£V=1+I+C_+-    •    •  (2.2) 

RT  V       y2 


where  B   and  C  are   referred   to  as   the    second  and   third  virial    coefficients, 
respectively.      This   equation  was   proposed   in   1885   by  Thiesen  to   describe    the 
behavior   of   gases   and  was    shown  by  Kammerl ingh-Onnes   (1901)    to  represent    the 
p-V-T  properties   of  many   gases  well.      Although   originally   empirical    in  origin, 
Mayer    (1937)    showed  that    this  expression  arose    naturally   from   a    statistical 
mechanical   model.      He  was  able    to  relate    the   experimental   values   of    the  virial 
coefficients   to   the  molecular  properties  of    the  material.      Typically,    the 
virial    equation  is   truncated  at   the    second  virial    coefficient;    such  an  expres- 
sion is   applicable   at    subcritical    temperatures  for  pressures  below   about   1   MPa 
and  for   increasing  pressures  at  higher   temperatures,    but   only  for   gases.      A 
major   advantage    of    the  virial    equation  of    state    is   that    second  virial    coeffi- 
cients are   known  for   a  huge    collection  of  materials    (Dymond  and   Smith,    1969) 
and  can  be   estimated   in  the   absence    of  measured  values    (Hirschf elder,    et   al.  , 
1954). 

The   reduced  equations  of    state   are   a    second  approach.      An  example    is   the 
Pitzer    (1957)    accentric   factor  method: 

Jj  =  Z°(Tr,Vr)    +  (o  Z1(Tr,Vr)  (2.3) 

The  functions  Z°  and  Z   are  universal  functions  of  the  reduced  temperature  and 
volume  derived  by  correlating  the  measured  properties  of  many  materials;  the 
accentric  factor,  w,  is  an  empirical  device  to  compensate  for  the  non- 


spherical    nature    of    the  molecule.      Conformal    solution  theory   (Henderson  and 
Leonard,    1971)    and   the  TRAPP  program  recently   developed  at  NBS   (Ely  and 
Eanley,    1981)    also  fall    into  this   category.      A  similar    idea   is   the    scaled 
equation  of    state   around   the   critical    point    (Stanley,    1971).      These    schemes   are 
based  on  the    idea   that  by    taking  molecular    size,    shape,    and  attractions   into 
account,    all  materials   can  be   compared   in  the    same  way.      These   forms,    although 
rooted   in  a    theoretical   model,    depend  on  data   for   real   materials   to   define    the 
functional   relationships.      Once    these   relationships  have   been  defined  by  a    set 
of  well-characterized  reference    fluids,    only   a  minimal    set   of   properties  are 
needed   to  estimate    the    scaled  properties   of   other   fluids.      Such  equations   of 
state   are   appropriate   at  all   conditions,    limited  only  by   the  knowledge    of    the 
reference   fluids  and   the   ability   of    the    scheme   to   compensate   for   the  various 
characteristics   of   real    fluids. 

One    of   the   largest   classes  of   equations  of    state,    comprising  a   large   fraction 
of    the   industrially  used  expressions,    are    the  modified  van  der  Waals   equations. 
The   original    equation  was   proposed   in  1873  by  van  der  Waals: 


'-.^t:.£ 


The   first    term   compensates   for    the   excluded  volume    in  a  bimolecular   collision; 
in  the    classical    interpretation,    b   is  four   times   the  molecular  volume.      The 
second  term   accounts   for    the   attraction  between  molecules.      Other   equations    in 
this    class,    such  as   the  Relich-Kwong-Soave    (1972)    equation,    and   the   Peng- 
Robinson   (1976)    equation  include   modifications  of    the   second  term,    and  use    it 
to   compensate   for   deficiencies   in  the   first   term   in  addition  to  accounting  for 


intermolecular   attraction.      These   equations  are   cubic   in  volume,    which   is   the 
lowest   order   of   equation  able   to  represent  both  liquid  and  vapor   behavior  with 
the   same  function.      Although   this   class  of    equations  of    state  has   a   liquid 
branch,    the  molecular  volume  has  been  derived  for   a   gas-like   fluid  and   thus 
there   is  no  reason  to  expect    that   these   equations  will    accurately  produce   the 
dense    fluid   (e.g.,    liquid)    conditions.       Indeed,     Henderson   (1979)    has    shown 
that    such   expressions  are  fundamentally   flawed  at  high  densities  because   the 
packing  problem  has  i.ot  been  properly  addressed.      Using  an  expression 
involving  a   large    number   of    arbitrary  parameters  can  improve  the  fit   in  the 
dense   fluid  region  but  has   two   disadvantages:      first,    it   gives  rise   to   calcu- 
lation complexities  and  second,    the  physical   meaning  of   the  parameters  is 
lost. 

The  final    approach   consists  of    equations  of    state   based  on  a    theoretical    fluid 
which  is  used  as  a  reference  to  which  real    fluid  behavior   is    compared.      In  a 
sense    this  approach   is   the   theoretical    analog   to  the  empirical   reduced 
equation  of   state.      The  most   important  feature    of   a  model    fluid  is   not   that   it 
describe   a   particular  real    fluid  exactly  but    that   it   contain  the  proper    treat- 
ment  of    the  high  density   states.      The  Leonnard-Jones   fluid  (Verlet   and 
Levesque,    1967)    is   an  example;    it   is   composed  of   molecules   that  repel    one 
another  at   short   distances  and  attract   one   another   at  moderate  and  large 
distances.       A  second  class   of    reference   fluids   are   those    characterized  by   a 
hard-body  repulsion.      The    state    of    development    of    this    class   of   equations   is 
demonstrated  by    the   fact    that   the  most  recent   international    steam   tables 
(Haar,    et   al. »    1984)    are  based  on  a  hard    convex-body    equation   of    state. 


A   simple,    well-characterized  hard-body    reference   fluid,    and  the   one   on  which 
the    equation   of    state   described  in   this  report   is   based,    is    the  hard-sphere 
fluid.       Such   a   fluid  exhibits   an  infinite   repulsion  force    for    a  bimolecular 
collision  at    some    distance    of    closest   approach;    there   is  zero   net    inter- 
molecular   force    at   greater   distances.      There  have  been  extensive   computer 
simulations   of    this   fluid   (Barker  and  Henderson,     1971)   and    statistical 
mechanical    treatments  to  describe    it    (Reiss,     1965).      There   are    three   analy- 
tical   functions   associated  with   the  hard   sphere    fluid:      the   Percus-Yevisk 
pressure   equation  (Lebowitz,    1964),    the  Percus-Yev  ick  compressibility   equation 
(Verlet,     1964),    and   the    Carnahan-Starling    (1969)    pressure    equation.       The 
Carnahan-Starl  ing  form    shall   be    used   in   this  discussion  primarily  because   of 
its  familiarity  to  the  engineering   community: 


£V  „    1   +  y  ±  v2   -   y3  (2  5) 

RT  (1  -  y)3 


where  y  =   b/4V. 

Although   the   Carnahan-Starling  equation  more   closely   represents  the  properties 
of    the  hard   sphere   fluid  at  highly   compressed,    nearly   close-packed   states,    all 
the  hard   sphere   fluid  representations   are   equally  good  at  lower,    more 
physically  realistic    densities. 

This   equation  was   derived  as   a   closed  form   solution  to  an  infinite   geometric 
series  which  was   in  turn  based  upon   the  virial   expansion  of  Ree   and  Hoover 
(1964).       Carnahan   and    Starling    demonstrated   that   eq.     (2.5)    accurately 
describes   rigid   sphere    behavior.       In    this    equation,    'b'    is    a   molecular  volume. 
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Although   it   has   a    similar    interpretation  to   the   b   in   the  van  der  Waals   or 
Redlick-Kwong-Soave    equations,    it  will    not  have    the    same   numerical   value.      In 
the   limit    of   large   volumes   this   expression  leads   to  a  van  der  Waal  s-1  ike    term. 
The    converse,    however,    is   not   true;    at    small   volumes   the  van   der  Waals 
equation  does  not   lead   to  eq.    (2.5). 

The  hard-sphere   equation   is  modified   to   compensate   for   long-range    attractive 
forces  by    the   addition  of   a    second,    semi-empirical    term  to   arrive   at   the 
equation  of    state    treated   in  this  work: 


pV  _    1  +  v  +  v2  -   v3   _ 
RT  '  ,,   _      v3 


(1  -   y)3  '  RT(V  +  b) 


(2.6) 


This   form,    first    suggested  by  DeSantis,    et   al.    (1976),    was   proposed  for    the 
prediction  of  multi-component,    vapor-liquid   equilibria   and  exhibited   good 
agreement   for    simple    inorganic  molecules,    hydrocarbons,    and   their  mixtures. 
(The    procedures   for   extending   this   equation  to  mixtures   are    discussed   in   the 
next    section.) 

The   physical    significance    of    the   parameters   defined  by    the   Carnahan-Starl ing- 
DeSantis   equation  of    state   and   their    temperature    dependence    can  be   understood 
by   comparing   the   differences  between  the   reference    and  real    fluids.      The   real 
fluid   is  not    composed  of  molecules  with   a  hard-sphere   repulsive    core.      The 
determined  value    of    the    'b'    parameter    is   an  effective    core    that    is   a  measure 
of    the    average    closest   approach   of  molecules.      In  the   real    fluid,    as   the 


temperature    is  raised,    the   average   kinetic   energy   also   increases  and  the 
average   distance   of    closest  approach  becomes   smaller*,  thus,    the  molecular    size 
parameter  will   decrease  with  increasing  temperature  because  of  the   softness  of 
the    short-range  repulsion  in  the    intermolecular   potential.       The    parameter   'a' 
has   the  dimensions  of   energy  per  mole   times  volume   per  mole*,   in  a   spherical, 
non-polar  molecule   it  would  be   connected  to  the   dispersion  forces.     The 
attractive   part   of    the  real    intermolecular  potential   can  be  highly  directional 
because   of   dipole  moments  and  higher-order   internal   charge   distributions.     As 
the   temperature    is  raised,    the    strong  coupling   of    these   directional    forces   is 
'washed  out'   and   the   average   attractive   force    drops.       The    temperature    depen- 
dence   of    'a'    is   a   reflection  of    the   non-spherical    character   of    the   force 
between  molecules.       In  this  work   the    temperature    dependence    of   'a'   and  'b'    are 
represented  by    the  following  forms: 

a   =    a^xpUjT  +    a2T2)  (2.7) 

b  =  bQ   +  bxT  +  bjT2  (2.8) 

These   forms   are   those    suggested  by  DeSantis,    et   al.    (1976)    with   the   addition 
of    the   T2    terms.      The   methods   of    determining   'a'   and   'b'   as   functions   of 
temperature   from    experimental    p-V-T  data   are   discussed   in  Section  6. 

This  hard-sphere   equation  of    state  was   selected  for    the   present    study   of 
refrigerant   properties  because    it  has   the   advantages   of   being  a   simple  expres- 
sion based  on  a   good  physical   model.      With  a  few   exceptions,    the    species   used 
as   refrigerants  are    small,    nearly   spherical    molecules.      These   attributes   in 
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turn  allow    only   two  physically   significant  parameters  to  represent  both  the 
liquid  and  vapor  phases.      The   equation  defines   a  fundamental    thermodynamic 
relationship  which,    combined  with  perfect   gas  heat  capacity   information, 
allows  a   complete   set   of   consistent   thermodynamic   properties   to  be   evaluated. 
The    equation  of    state   can  also  be   extended  to  mixtures. 
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3.      THE  EQUATION  OF   STATE  AND  MIXTURES 

In   the   previous   section,    we  have   discussed  the   shortcomings   of    equations   of 
state  with   the  van  der  Waals  excluded  volume    term  and  how    the   equation   of 
state   for    the  hard   sphere  fluid  rectifies   this   serious   physical    flaw   of    the 
van  der  Waals  model.      In   this    section,   we    shall    discuss  how    the   equation  of 
state   for  mixtures   is  produced. 

The  first   approximation  made    in  the  application  of    the   CSD  equation  of    state 
to  mixtures   is   the   assumption  that   no   substantive    change   needs   to  be  made   to 
represent   the  properties  of    a  mixture.      That  is  to   say    that   there   exists  an 
effective    'a'    and   'b'    such    that   an  equation  of   form  identical   to   the   one   used 
for   pure  materials   can  be    used  for    the  mixture.      This   section  will    focus 
primarily  on  the   evaluation  of      such  an  effective    'a'    and   'b'.      A  different 
approach,    accounting  for    the  behavior   of   a  mixture   of   hard   spheres  of 
differing   diameters,   will   be  mentioned  briefly  at   the    close    of   this    section. 

In   the    simplest  mixture  models,    the   effective  molecular  parameters   are   defined 
as  follows: 


=    X     L    xixjaij  (3*1) 

i=l    j=l 


and 


-  E  £  *i*jt>ij  <3-2> 


i=i  j=i 
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When  i  =   j,    the  values   of    a.  .    and  b^   are   those    of    the   pure  materials.      The 
values   of   &z-    and  b::    can  be   obtained   if   nearly   any   experimental    property   of 
the  mixture    is  known;    examples  will   be    given  later    in  this   discussion.      The 
motivation  of   a  mixing  rule    scheme   is   to  determine    the  values   of   a^.    and  hz- 
without   the   assistance    of   measurements.      The   techniques   used  for   predicting 
the  mixture   parameters  fall    into  at  least   two   classes,    those  rooted   in  a 
rigorous   or    at   least   a    semi-rigorous   physical    model    and   those   rooted   in  the 
correlation  of  large    sets   of  mixture   data    (Fredinslund  et   al.    1977)*,    the 
latter    schemes   are    typically  not   applied   to   species  with   one    carbon. 


The    origin  of    the    size   parameter,    b12>    *s   tne  morc   obvious   of    the    two.      The 
closest   distance    of    approach  for    two  hard   spheres   is   the    sum   of    the    two  hard 
sphere   radii  which  is   evaluated  as   follows: 


biz  =  — — r* — -  (3-3) 


In  this  Note,    Eqn.    3.3    is  approximated  as: 

b12  =    (bx   +  b2)/2  (3.4) 

This  approximation  has  been  made    primarily   for  mathematical    simplicity.      When 
the  values  for    'b'    are    similar  for   the    two   components,    the    two  mixing  rule 
schemes   can  barely  be   distinguished.      Unless   the  molecular  volumes  differ  by 
more    than  42  %,    Eqn.   3.3    and  Eqn.   3.4  will   not   differ  by  more    than  1  %*,    and, 
hence,    will   not   give   rise    to  more   than  a  1/2  %  difference    in  the  value    of  b. 
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The   physical    origin  of    the  mixing   scheme   for   a,2    has   its  origin  in  a  variety 
of  molecular   and  quantum  mechanical   models.      In  the   equation  of    state,    the  value 
of    'a'    arises  from   the   attraction  between  the  molecules.      On  a  molecular 
level,    the   interaction  energy   is   typically  expressed  as  a   correction  to  the 
geometric  mean  of    the   pure   component  energy   force    constants: 


"12--   (1  "   f12)(alla22)1/2  (3'5) 


The   correction  factor,    (1  -    f-2)  »    attempts   to  account   for   differences   in 
polariz ability  of    the    two   species    (the   ease  with  which   the    two  electron  clouds 
are   deformed)    and   the   closest   distance    of    approach.      Schemes  for   predicting 
f-,    are  reviewed  by   Pesuit    (1978),   who   considers   them  on   several   levels   of 
empiricism.      What   is   impressive    is   the   degree    to  which   the    interaction 
parameter   can  be   estimated  for  relatively   complex  molecules.      Equally   impres- 
sive however    is   the   unexpected  collapse    of   these   estimation  rules.      Enobler 
(1978)    has    shown  that   even  for   so    simple  mixtures  as  He/Ne,    He/Ar,    He/Kr  and 
He/Xe    that   an  estimation  of   the   interaction  parameter   such   as   the   following   is 
far  from  adequate: 


h   +   h        (<7U  +  a22)6 


I.    is   the   ionization  potential    of   component    i;    and  a ..,    the  molecular 
diameter.      The   experimental  value   for    (1  -   f,2)    *s  0.611   for  He/Xe;    its   pre- 
dicted value    is  0.824.      The   conclusion  that   one  would  draw   from   the   examples 
given  by  Enobler   and  Pesuit  are    that   the  methods  for  estimating   fj~    are    not 
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very   satisfactory   and   that   a  method  appropriate   for  one   class  of   mixtures  may 
not   be    satisfactory   for    another. 

The    thesis  of    these   models    is   that   all    the    information  needed  to  describe   the 
interaction  between  unlike   species  can  be   derived  from    information  about   the 
interaction  between  like   molecules.     This   thesis  neglects   the   possibility   that 
there   can  be  interactions  in  mixtures   that   are   either   unimportant    or   do   not 
occur    in   the   pure   fluids;    the   converse    is  equally  possible.      Recently,    Gubbins 
and  his   colleagues  have  attempted   to  evaluate   functions  related   to   fj~    from 
the  mul  t  i  pole- mult  i  pole  interactions  between  molecules  (Gubbins  and  Twu, 
1978).      Wallis  et   al.    (1984)    have    shown  that    this  approach   can  be    used   effec- 
tively   in  describing  mixtures  of   carbon  dioxide   and  ethane.      In  general,    one 
should  approach  these  estimation  schemes  with  caution  and   in  the    spirit   that 
they  can  be   used  as   guides  when  no  other   information  exists.      One    should  be 
aware    that,   at   present,    the   only  quantitative  estimation  of   a   function   such  as 
fj~    c&n  arise   from  experimental    information. 

In  this  Note,    the   interaction  parameter  will   be   evaluated   solely  from 
experimental    mixture   information.      Almost  any   small    set  of  data   on  the  mixture 
is   sufficient   to  evaluate   f -. 2*      Because    there    is  but  a    single   parameter  to 
evaluate   for    the  mixture,    one   needs    in  principle   only   a   single  measurement; 
however,    the  larger    the   experimental    data    set,    the    greater   the    certitude    in 
the   evaluation  of   f.,.       In   section  6   of    this  Note,    a  description  of   the   use    of 
bubble   point   information  to  evaluate    f.2    *-s   given.      One    often  finds   cross 
second  virial    coefficients,    B,2»    in  the   literature;    such   information  can  be 
used  to  find   iyy    through   the   appropriate   expansion  of    the    equation   of    state. 
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One    can  also   use    gas   phase    data.      Of    the   available   data,    liquid  phase    densi- 
ties are   perhaps    the   least  reliable   because   a    tiny   difference   between  the 
actual    density   and   the   equation  of    state   density    at    the    same   pressure   can  lead 
to  large    errors    in   the    pressure. 

The    above    discussion  has    dealt  with   finding   an   effective    'a'    and   'b'    for    a 
mixture   for   use   with  an  expression   developed  for   pure    fluids.      There    are 
equations   of    state    that   approximate   the   composition  dependence   for  mixtures  of 
hard   spheres    (Mansoori   and  Leland,     1970)    such  as   the   following: 


pV  =    1  +   (1  -  3A)y  +   (1  -  3B)?2  -   (1  -   Ov3 
RT  (1  _   y)3 


(3.7) 


In   this   equation,    the  various   terms    in   the   Carnahan-Starl  ing   equation  have 
become   functions  of   the   averages   of  various   powers   of    the   hard   sphere 
diameters.      When  the    two  molecular  volumes   differ  by   a   factor    of   1.5,    A  is  no 
larger  than  0.018,   B  no  larger  than  .0089,   and  C  no  larger  than  0.013.      In  the 
limit   of    a   pure   material,    A,    B,    and  C,    become  zero  and  Eqn.    3.7    reduces   to   the 
original    Carnahan-Starl  ing  expression.      In   this  work,   we    shall    assume    that   the 
departure   from    the   Carnahan-Starl ing   equation  of    state    caused  by   the  disparity 
between  the  molecular   sizes  is  unimportant  and   that    this   disparity   can  be 
compensated  for    in  the    interaction  parameter   f.. «. 
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4.      THERMODYNAMIC  FUNCTIONS  AND  REFERENCE   STATES 

There   are    two   themes  found   in  each   of    the    two   parts   of    this    section.      The 
first    is  a   general    description  of    some   important   thermodynamic   derivatives. 
In  this  Section  on  pure   fluids,    the    derivation  of    the    thermodynamic   properties 
from   a  p-V-T  equation  of    state    is   discussed.      In  the   part   on  mixtures,    the 
ideal   mixture  model    is   developed.      The    second   theme    confronts   the    problem   of 
reference  materials  and  conditions    (or   reference    states).      In  the   section  on 
pure  materials,    the   perfect    gas  appears  as  a   necessary  part    of    the    derivation 
and,    through   it,    other   reference    states  are   constructed.      In   the    section  on 
mixtures,    the  reference    states   of    the    components   of    the  mixture   are    connected 
to  the   thermodynamic  properties  of    the  mixture. 

Pure  Fluids 

In  the  previous   sections,    we  have   discussed  the  physical    origin  of    the  hard 
sphere   equation  of    state  and  the  merits   of   using   such   a  model   at   high  fluid 
densities   instead  of   a  model   where   the   excluded  volume    is   expressed  by    the 
van  der  Waals,    l/(V-b),    term.      We   begin  here   a  review   of    the    derivation  of    the 
thermodynamic  properties  of    a  fluid  from   its  p-V-T  equation  of    state   and   show 
the   specific  results  for   the    Carnahan-Starling-DeSantis   equation.      In  this 
discussion  we  will    find  that   the   pressure   equation  of    state   does  not   include 
all    the   information  necessary   to   evaluate    the    thermodynamic   properties   com- 
pletely.     We    shall    resort   to   the  properties  of    the  perfect   gas,    spectroscopic 
information,    and  relationships   provided  by    statistical    thermodynamics   to  form 
a   complete    set   of    thermodynamic   functions. 

The   pressure   is    connected  to   the  Eelmholtz   free   energy  and,    through   it,    to   all 
the   other   thermodynamic   properties  by    the   following  relation: 
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»  -  -ft)  T  «•" 


The  Helmholtz    free   energy  can  be   evaluated  from    the  pressure  by   an  integration 
over  volume.      Several    difficulties  are   encountered  in  this   integration:      first 
the    choice    of    appropriate   limits  for    the   integration*,  second,    the   possibility 
of   divergences   in   the   integral*,  and,    finally,    the   recognition   that    there  will 
be    integration  constants.      The  problems   of    integration  limits  and  divergences 
are    closely   connected*,    they  will   be   discussed  together. 

When  the  Helmholtz    free   energy   is  calculated,    we  would  like    to  choose    the 
limits  of    the   integration  so   that  one    of    them   is   independent   of    the   form   of 
the    equation  of    state.      This   can  be   achieved  by   having  one   limit  be  V  =   eo, 
where   all    gases  become   perfect    gases*.      At   that  limit,    however,    the    integral 
diverges.      This   problem   can  be  resolved  by   evaluating   the    difference   between 
the   properties   of    the   fluid  described  by   the   equation  of    state   and  those    of 
the    perfect    gas.      As   the  volume   approaches   infinity,    this   difference    becomes 
zero    sufficiently   strongly   that   the   integral    remains   finite.      Thus, 

P"P°     -  P  -  f  -  -(L(A  -  A"))  (4.2, 


*The    term   'perfect    gas'    is   used   rather    than  'ideal    gas'    in  order    to   preserve 
the    term   'ideal'    for    certain  properties   of  mixtures. 
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and 


A(V,T)   -  A°   (V,T)    =  /   (p  -  ^)dV  (4.3) 

V 


Throughout    this    section,    functions  applicable    to   the   Carnahan-Starling- 
DeSantis  equation  of    state    (such   as   those    evaluated  by  Equations  4.2    and  4.3    ) 
are    given  in  Table  4.1. 

Before   evaluating  the   remaining  functions,    let   us   review    the   features   of    the 
Helmhol  tz   free   energy   calculation.      First,    in  performing   the   integration,   we 
have   resorted  to  a   reference   material,    the   perfect    gas.      It  was   chosen  to 
avoid  a   divergent   integral.      Although   the   perfect    gas   is  a   construct,    it   is  a 
convenient   and  useful    material    against  which   to   compare    the   properties  of    any 
other  material — all   gases  become   perfect   gases   in  the  V  =  °°  limit   and   its 
properties   are   known  to   the  last   detail.      The    second  feature    of    the    integra- 
tion  is  a  more   general    consequence    of    that   operation.      Since    the    integral   was 
over  volume,    both  A(V,T)    and  A°(V,T)    may  be   in  part   functions   that   are    inde- 
pendent  of  volume.      Because    the   difference   between  the   two  free   energies    is 
uniquely  defined,    the    temperature   dependent  and  volume   independent   parts   of 
A(V, T)    and  A°(V, T)    must   be    identical.      Unless   one    is  aware   of    these   functions 
which  are    constants  with  respect   to   the  volume    integration,    other    thermody- 
namic functions  will    not   be    complete. 

Once    the  Helmhol  tz    free   energy  has  been  evaluated,    all    the   other   thermodynamic 
functions  arise   by   the   usual   operations 
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5  =  -  {  2A)  (4.4a) 


S-    S°=  -(L     (a-  A0))  <4'4b> 

\  3T  /y 


and 


E  =  A  +  TS  (4.5a) 

E  -  E°  =    (A  -  A°)    +  T(S  -   S°)  (4.5b) 


and 


H  -  E  +   pV  (4.6a) 


H  -   H°  =    (E  -   E°)    +   (pV  -  RT)  (4.6b) 


and 


G  =   A  +   pV  (4.7a) 

G  -  G°  =    (A  -  A°)    +   (pV  -  RT)  (4.7b) 

The  values   of    the   functions    in  Equations  4.4b,    4.5b,   4.6b,    and  4.7b   are    given 
in  Table  4.1.      The   reference  material    for    the   functions  enumerated   in  these 
equations  has,    in  all    cases,    been  the   perfect    gas   at    the    same    temperature   and 
volume   as   the   real    fluid.      For  E  and  H,    the    perfect    gas  values   depend  only   on 
temperature,    hence,    volume    is   deleted  as   an  appropriate    argument   for    those 
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properties    in  Table  4.1.      The   reference   condition  for   the  Gibbs   free   energy    is 
often  chosen  terms   of   a  fixed  pressure.      This   choice   arises  because   of   the 
role   of    the  Gibbs   free   energy    in  determining  phase   equilibrium.      Equation  4.7b 
can  be   easily  modified   to  account  for   a  reference   pressure,    p*,    as  follows: 


G(V,T)   -  G°(V,T)    =  G(V,T)    -  G°/p  =  ^,    T 


=  G(V,T)   -  G°(p,T)    +  G°(p*,T)    -  G°(p*,T) 


=  G(V,T)    -  G°(p*,T)   -  RT  ln^ 


1  p*y 


or 


G(V,T)    -  G°(p*,T)    =  RT  £n  (%&  )  +     G(V,T)   -  G°(V,T)  (4*8) 

\P*V/ 


This  function  can  be   found   in  Table  4.1. 


The  heat   capacity    at   constant  volume    is   evaluated  as  follows: 


cv  =  (af)v  =  T(fi)v  (4-9a> 


C     -   C°  =  (Z-(B  -   E°))      =   t($-(S  -    S°))  (4.9b) 


Earlier   in   this    section,    the    need   to  recognize    the   existence   of  volume 
independent   functions  as   a   consequence    of    the   calculation  of    the  Helmholtz 
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free   energy   from    the  pressure  was  noted.      It   is   in  the  evaluation  of    the  heat 
capacities   that   the   existence    of    these  functions   becomes  most  apparent.      The 
perfect   gas  has  both   internal    energy   and  a   heat   capacity,    yet  both   are 
independent    of    the  volume. 


(dE°)     =  T(dj>°)  _   po  =  RT  -   RI  =  o 

Vav  /T       VaT  A,  v       v 


(4.9a) 


/dc°\   _    aV     _    aV     .  T(d2v°\     =0  (4.9b) 

\3V  /T  3V3T  3TdV  \dJ2   )y 


The   perfect   gas  values   of  E  and   Cy  as  well   as  H  and  C     and  parts   of   A,    G,    and 
S   can  be   evaluated  from   calorimetric  measurements*,  however,    these  quantities 
are    typically   calculated  through   straight-forward   statistical   mechanical 
relationships  from    infrared  and  Raman  Spectra   of    low   pressure   gases. 

The  heat   capacity  at    constant   pressure    can  be    calculated  from    the   enthalpy   or 
entropy   at   constant   pressure 


c>  ■  (§?)  -  T  (& 


-/p  (4.10a) 


however,    because   the  equation  of    state   is   a  function  explicit   in  temperature 
and  volume,    C     is  more  readily  evaluated  from  Cy   and  temperature  and  volume 
derivatives   of    the    pressure: 
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s  ■  °v  - T  (E)V/(K)T  <"*> 

From   the  beginning  of    this  derivation  we  have  been  forced  to  choose   a 
reference  material.      The   initial    choice,    the   perfect    gas,   was  made    to   avoid  a 
divergence    in  the  value   of   an  integral.      There  may  be   reasons  for   choosing 
other  reference    conditions.      Some    properties  require   no  reference,    properties 
such  as  pressure,    volume,    temperature,    and  heat   capacity  can  be  measured 
absolutely.      Indeed,    the  reason  for  using  the   perfect    gas  reference  for  heat 
capacity    is   that   it   can  be   determined  quite   accurately.      Other   properties  have 
only  relative  measures;    a  zero  value  has  a    significance   only  by   convention. 
For    these    properties,    only   changes   can  have   any   physical    significance.      Among 
these   are   internal   energy,    enthalpy,    and  entropy.      The    two  free  energies 
present   a  very   special    problem.      A  change    of  G  when  temperature    is   kept 
constant  requires   knowing   the  volume.      A  change    of  A  at   constant   temperature 
requires  the  pressure,    obviously   there   is  no  problem   in  either   case.      When  the 
temperature   is   changed,   however,    the  value   of   S  must   be   known.      Depending  upon 
ones   interpretation  of    the    statement   of    the   third  law   of    thermodynamics,    one 
can  debate  whether    the   entropy  has  a   known  absolute  value    or   not,    surely   its 
value   relative   to  a    special    state   at  0  K   is  known.      Thus   changes   in  G   and  A 
cannot   be   calculated  without   carrying  an  unknown  quantity  T  S(0  K)  . 
Fortunately,    one   uses  G   and  A  only    in  fixed  temperature   processes.      We   have 
then  at  least   three    kinds   of    thermodynamic   functions:      those   known  absolutely, 
those    known  only   relative    to  a    state   at   the    same    temperature;    and   those   known 
relative   to  a    state   at  a   different   temperature. 
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One   may   have    special    reasons  for    choosing  a   particular   reference    condition. 
One  might   choose   a  readily  accessible   physical    state.      For   example,    the   refer- 
ence   state    in  chemical    thermodynamics    is   often   the   element    in   its   stable   form 
as  found  at  101.33   kPa   and  25°C.      In  electrochemistry,    all   voltages   are 
typically   compared   to   the    standard  hydrogen  electrode.      For   gases,    one  might 
choose    the   perfect    gas   at  101.33   rPa   and  25° C.      For   refrigerants   the   reference 
is   typically   the   liquid  at    its   saturation  pressure   at  -40°C.      There   are,    of 
course,    as  many  reference    states  as   there   are  applications  for   thermodynamics', 
they    are   all    interconnected   through   thermodynamic   relationships.      Some   of 
these   reference    states   carry    the  weight    of   an  internationally   accepted 
definition  and  are   often  referred  to  as  STANDARD  STATES.      The    significance    of 
a  STANDARD  STATE  is   its  widespread  use   and  precise    definition. 

The   next   few   paragraphs  will   be    devoted  to  examples   in  which   the   enthalpy, 
internal   energy,    and   the   entropy   are   calculated  with  respect   to  an   arbitrary 
reference    temperature,    T     f,    and  reference    pressure,    Pref.      Since    the   pressure 
is   an  explicit   function  of  volume   and   since    several    states   of   a  material  may 
have    the    same   pressure,    the   description  of    the   reference    state    is  further 
defined  by   indicating   the  volume    parenthetically   after    the    pressure. 

The  most   effective  way   to   go   from  one    thermodynamic    state    to  another    is   to 
break,   the    path   into    steps  for  which  a    single  variable    is    changed  at  a    time   and 
for  which   the   changes   are   readily  calculated.      For   example,    for    the   enthalpy 
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H(p(V),T)   =  H(p(V),T)   -  ff>(T) 
+  ff»(T)   -  H»(Tref) 

+  H°<Tref>    -H<Pref<Vref)'   W 

+H<Pref<Vref>'    Tref>  <4;"> 

The    terms   in  the  first  and  third  lines  on  the   right  hand   side   of  Equation  4.11 

can  be   evaluated  according  to  the  enthalpy   expressions   in  Table  4.1.      The 

o 
second  line    involves   the   integral    of   C     between  the   two  temperatures,    T  and 

T__«*      Th©  final    term  is   the  reference  enthalpy,   which  may  be    set   to  zero  by 
re  i 

convention  or    is  more   properly  moved  to  the  left  hand  side    of    the   equation. 
We  have    then 

H(p(V),T)   -  H(pref<Vref),    Tref)    = 
[H(p(V),T)   -   H»(T)] 

T 

r     o 
+  /  Cp  dT 

Tref 
"    [H(Pref(Vref>'   Tref>    "  H°<Tref)]  «.12> 

Internal   energy  can  be   evaluated  in  a    similar  fashion 

E(V,T)   -  E(Vref,    Tref)    = 
[E(V,T)    -  E»(T)] 

T 

r     o 
+   7    Cy    dt 

Tref 
~[E(Vref,    Tref)    -  E«(Tref)]  (4.13) 
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As   in  the   case    of   enthalpy,    the  first  and  third  lines  on  the  right  hand  side 
of    this  expression  are  found  in  Table  4.1. 

The   derivation  of   the  entropy    is   similar,    however    it  has  an  additional    term 
because,    unlike  H°    and  E°,    S°    is  not   independent   of  volume. 

S(V,T)   =   S(V,T)   -   S°(V,T) 

+  S°(V,T)   -   S»(Vref,T) 

+  So<Vref,T)   -   So(Vref,    Tref) 

+  S°<Vref    Tref>   "    S<Vref    Tref> 

+   S(Vref    Tref>  <4'14> 

In  this  expression,    the  first  and  fourth  terms  are    gr.ven  in  Table  4.1,    the 
third  term   is  the   integral    of   c£(T)/T,    and  Vfe£    is  the  volume  at   the  reference 
state    (e.g.,    for  refrigerants,    the    saturation  volume   of    the  liquid  at  Tre£> . 
The    second  term    is   evaluated  as  follows: 

oo  Y       /dS°\ 

S   (V,T)   -    S  (Vref,T)    =  /       [Jy-)      dV  (4.15) 

V     *  T 

Tref 


The  Maxwell    relation  allows   this   derivative,    (dS/dV)m  to  be   evaluated  in  a 
straightforward  manner. 
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(-)  =  -  0 


3_b\  (4.16) 


f  n3T  V 


For    the   perfect    gas: 


\3TA,  V 


The    second  term   then  becomes 


S»(V,T)   -   S«(Vref,T)   =  -  /       (^)dV  =-Rln    J 

Vref  V    ref 


The  quantity    that  appears   in  the   refrigerant   property   tables   is  as  follows: 


S(V,T)   -   S(Vref,    Tref)    = 


[S(V,T)   -   S»(V,T)]   -   R  ln(  ^ 


Vref 


T 

+       /       (CV  )  dT 


Tref V  T 


"    <S<Vref    Tref>    "   S°<Vref>   -   Tref)]  (4.17) 
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Table  4.1      Thermodynamic  Functions  Arising  from    the  Perturbed  Carnahan- 
Starling  Equation  of   State 


pV  _   1  +  y  +  y2  -  y3 a 

RT  (1  _  y)3  RT(V  +  b) 

y  =  b/4V 

A(V.T)    -   A°(V,    T)    =-iln  (£+-*-)   +     4RTP        +  — *^ 

b         \      V     J        (V-   B>         (y_  f 

p  =  b/4 

G(T,p(V))    -  G°(p*,T)   =  RT  m  (*£-)  -  ■»  I  n(W  +  b)  +  — ^ (8V2  -   9Vp   +  3p2) *- 

Vp*V/        b       I      V    )        (y  _  p)3  H  V  + 


(V  -   B)2 


lid.pOD.^)    =  GCT.pm.Xi)    +   (1  -   ii)(|^ 

Xi   T.p 


=  ^(,..I.,4)    +  RTUW  +   RTP(4V  -  3P)    +  RI^j(4V2  -  2VP)    +  ab^        /v^ 
1  W  (V     -   P)2  (V     -  6)3  b2  V     V 


abi         +   2*iai  +  2*iaii   £     /     V 
b(V  +  b)  b  \V  +  b 


Note:      the  above   expression  for   chemical    potential    applies  only   to  a   binary  mixture. 

S(V  T)   -   S°(V  T)    =   a'b  ~   ab'   I  JW  ±  b\  +  ab'        -  RB(4V  -  3  6)    _  RTB'UV2  -  2VB) 

b2  V     V     y        b(V  +  b)  (y_  p)2  (v_   p)3 

E(V  T)    -  E°(T)    =    a'bT  "   ab'T  "   ab  :  n(W  +  b) i         ab'T       _   RT2B'(4V2  -  2VB) 

b2  V     V     /      b(V  +  b)  (v  _   p)3 

H(T  p(V))   -  H°(T)    =    ab'T  -   a'bT  -   ab  „  p/V  +  b\  ,    ab'T-   ab  +  RT(4V2  -  2Vp?(p  -   p'T) 

b2  V     b     )     b(V  +  b)  (y  _   p)3 

Cv(V,T)   -  c£(T)    =  gRI^B^CVB  -  2V2)    +  2RTV((B"T  +  2B')(B  -  2V)    +   B'2T)    _        Tab' 

(V  -  p)4  (V  -  P)3  b(V  +  b)2 

+  T(ab"b  +  2a'b'b  +  2ab'2)    _    ( a"b2T  -  2a'b'bT  -i-  2ab'2T  -    ab"bT)   f  ^V  ±  b 
b2(V  +  b) 


cp=cv-t(S)2     /(Si) 


V. x  T. x 

NOTE:      primed  quantities   indicate   a    temperature   derivative;    double   primed 
quantities  are    second  derivatives  with  respect   to   temperature. 
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The   Thermodynamic   Properties   of  Mixtures 

The   above    discussion  presented   in  detail   for   a   pure    substance:       1)    the 
calculation  of    thermodynamic   properties  from    an  equation  of    state,    2)    the 
perfect    gas    state,    3)    reference    states. 

In   this    section,     the   discussion  will  be   expanded  to  include    the   thermodynamic 
properties  of  mixtures.      The   discussion  of  mixtures   in  many  ways  will    parallel 
that  for   pure  materials  but  will   differ   in  two   important  ways  because   of    the 
extra   degree(s)    of   freedom   accessible    to  mixtures,    the    compositions.      The   first 
topic    to  be    discussed  will   be   reference    states  for   mixtures.      Integral    to   that 
discussion  w  ill   be    the   description   of   a   special    class   of  mixture,    the    ideal 
mixture.      The    second   topic   that  will  be   covered  is  the  class  of   thermodynamic 
properties,    the   partial   molar   properties,    that   arise   from    the    derivatives   of 
extensive   properties  with   respect    to   changes   in  the   composition.      In   this 
class  of   properties,    the   chemical    potential  will    be    needed  to   evaluate    condi- 
tions   of    phase    equilibrium. 

If    the   components  of    a  mixture   did  not   tend   to   separate  during  a   phase    transi- 
tion,   one   could  define  reference   states  for  mixtures   in   the    same  way    that    such 
states   are   defined  for   pure  materials.      To   use   refrigerant  materials   as   an 
example,     one    could  choose    the  liquid  phase    at   its  bubble   pressure   at  -40° C. 
The   incipient  vapor  in  equilibrium  with    that  liquid   typically  has  a    different 
composition.       If    that  vapor   were   condensed  and  its  reference   point  were 
defined  as   the   liquid  at    its   own  bubble   pressure  at  -40°C,    one  would 
implicitly   have    changed   the  reference   point    of    the    parent  liquid   phase;    chaos 
would  follow!      The   following  paragraphs  discuss  a   solution  to  this  dilemma. 
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The   properties   of   pure  materials  were    initially   tied   to   the   properties  of    the 
perfect    gas.      Through   that   connection,    the   properties   of    the  material   at  any 
two  arbitrary   states  are   automatically   connected.      The   properties  of   mixtures 
will   be    connected  to   the   properties  of    their   constituent   components   through 
the    ideal    mixture.      In  the   discussion  of    the    ideal   mixture,    we    shall    find   that 
the   ideal  mixture   of   perfect    gases  will   be   particularly  useful. 

The   properties   of    the    ideal   mixture   can  be   described  in  terms  of    the 
following   simple   experiment:      the  mixing  of    the   pure   constituents   all   at   the 
same    temperature   and  pressure,    to  form   a  mixture  at   the    same    temperature   and 
pressure.      Regardless   of    the    composition,    temperature,    or  pressure,    one  would 
find  for   an   ideal    mixture   that  first,    the   final   volume    of    the  mixture  would 
equal    the    sum   of    the  volumes   of    the    unmixed   components.      That   is, 

V.CT.p.Cxi))    -2>iVi(T,p)    =  o  (4.18) 

i 

(Where  V^    is   the  molar  volume   of   component    i   or   the  mixture   and  i-    is  the  mole 
fraction  of    the   respective    component   in  the  mixture.) 

Secondly,    the  mixture  would  have   to  be   neither  heated  nor   cooled   to  bring   it 
to   the    same    temperature   as   the   unmixed  components*,    the   enthalpy  of   the  mixture 
would  be    the    same   as   the    sum   of    the   enthalpies  of    its   constituent   pure 
components: 

yT.p.lxj})   -L»iHi(T,p)    =     0  (4.19) 


30 


By  multiplying  Eqn.  4.18  by  the  pressure  and  then  subtracting  the  result  from 
Eqn.  4.19  we  have  a  similar  relationship  for  the  internal  energy,  E. 

Em(T,Vm(p),{xi})  -X^x.EidW^p))  =  0  (4.20) 

i 

When  the   temperature  derivatives  at   constant   pressure   are   taken  of  Eqns.   4.18 
and  4.19,    one  finds   similar  relations  for   the   coefficients  of    thermal 
expansion  (weighted  by   the  molar  volume)    and  the  heat  capacity   at   constant 
pressure. 

VT.p.Uil)   Op^T.MxiH-I^iViCT.p)   opi(T,p)   =  0  (4.21) 


(Where  a     is  the   coefficient   of   thermal    expansion  of    the   respective   component 


or  mixture,   I  (^)    .) 
V  VaTA 


Cpm(T'P(V)'{xi})   -ExiCpi(T'p(V))    -  0  <4-22> 

i 


(Where  C     is  the  molar  heat   capacity   at   constant   pressure   of   the   respective 
component    or  mixture.) 

When  the   pressure   derivative  at   constant   temperature    is   taken  of  Eqn.   4.18,    a 
linear  relation  for    the   isothermal    compressibilities   (weighted  by   the  molar 
volumes)    can  be   derived: 
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Vm(T,V.Ui})fiTm(T,p.Ui})   -£xiVi(T,p)PTi(T,p)   =0  (4.23) 


(Where  pj  is  the   isothermal   compressibility   of    the   respective   component   or 
mixture,   "  ff     (&)  .) 

v   \ap/T 


A  similar  relationship  can  be   derived  for   the  molar  heat   capacity   at   constant 
volume.      The    simple   linear  relationship  is  valid,   however,    only  under  very 

restricted  conditions.      By  using  the   following  relation  between  C     and  C    : 

v  p 

Cy  =   Cp  +  TV  aJ/PT  (4.24) 

the  first   term  on  the   right-hand- side   is  linear   in  composition  for  ideal 
mixtures   (Eqn.   4.22).      The   second   term  can  be  made   linear  with  respect   to  the 
composition  only    if  both  the  numerator  and  the   denominator   are   insensitive   to 
the   composition.      Except   for   a   truly  unusual    liquid  mixture,    the   only  material 
that    satisfies   such  a  requirement  is   the   perfect    gas.      The   perfect   gas   thus 
plays  a   central    role   in  connecting  the  properties  of   the  mixture   to  those    of 
its   constituent   components.      For   perfect    gas  mixtures. 


CVm(T,V(p),{xi})   "ExiCVi<T'Vi<P)>   =  °  (4'25) 


(Where  Cy   is   the  molar  heat   capacity   at   constant  volume   for   the  respective 
component   or  mixture.) 
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There    is  one    further   property   of    the   ideal   mixture,    the   entropy   of   mixing. 

Sm(T,p(V),{xi})    -£*iSi(T,p(V))    =-R2xiln(xi)  (4.26) 

i  i 

By   comparing  Eqns.   4.19   and  4.26,    the  Gibbs   free   energy  associated  with  mixing 
is  as  follows: 

Gm(T,p(V),{xi})    =   H^LpCY).^})   -   TSm(T,p(V),{xi}) 

or 

Gm(T,p(V),{xi})    -^XiG^T^pCV))    -  RT^XilnCx.)  (4.27) 

i  i 

The   definition  of    the   properties   of    the    ideal   mixture   allow    the   real    mixture 
properties  and  the  real   pure    component   properties   to  be    connected,  whether   it 
be    the   pure   components  at   the    same    temperature   and  pressure    or    the   pure 
components   under    some   other   conditions,    in  particular,    their   respective 
reference    states.      The   details   for   determining   the  mixture   equation  of    state 
have   been  discussed   in  Section  3    of    this  Technical   Note.      Once    the   form   of    the 
equation  of    state  has  been  established,    all    the   thermodynamic   properties  of 
the   fluid  mixture  with  reference    to   the   respective    properties   of   a   perfect 
gas  mixture  with   the    same   composition   can  be    evaluated  by   using   the   formalism 
described  earlier   in   this   Section.      Two  examples  follow:      the   evaluation   of 
the   enthalpy   of    the  mixture   first  with   respect    to   the   enthalpy   of    the   pure 
components   at   the    same    temperature   and  pressure   as   the  mixture   and   then  with 
respect   to   a   reference    state    of    the   pure  materials.      For    the    sake    of    simpl  i- 
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city,    a  binary  mixture  will  be    discussed;    this   example  may  be   expanded  to  any 
number   of   components. 

H^T.pW.X!)    =  H^pCV)  ,xx)    -  H^Xj) 

+  ^(T^p    -  XjHjCT)    -   (1  -  xx)^(T) 

+  x^Hjm   -   H1(T,p(V1))] 

+   (1  -   Xl)[^(T)   -   H2(T,p(V2))] 

+  x1H1(T,p(V1))    +   (1  -  x1)B2(T,p(V2))  (4.28) 

By   inspection,    one    can   see    that   this  expression  reduces  to  the   identity 
Hm(T,p(V)  ,xj)    =  H^TxpCV)  ,xj)  .      The   first,    third,    and  fourth  line    of    the 
right-hand- side    of   Eqn.   4.28    can  be   evaluated  using   the   functions   in  Table  4.1 
and   the   appropriate  values  for    the   equation  of    state   parameters,    a  and  b.      The 
second  line   equals  zero  because   a   perfect    gas  mixture   is   ideal.      The   final 
line    is  the  weighted   sum   of    the   enthalpies  of    the   pure   fluids   at   the    same 
temperature   and  pressure  as   the  mixture.      The    sum  of    the  first,    third  and 
fourth  lines   in   the   right-hand- side    of   Eqn.   4.28    is   the   enthalpy   of   mixing  of 
the    components    '1'    and    '2'    at  fixed  T  and  p.      A  simple   expansion  of    the 
calculation  that   led  to  Eqn.    4.28   produces   the  value    of    the   enthalpy   of    the 
mixture  with  respect   to   the  reference   pressure  and   temperature    of   its   consti- 
tuent  components,    as  follows: 
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Hm(T,p(V),Xl)    =   Hm(T,p(V),x1)    -   ^(T.Xj) 

+  H^er.xp  -  x^m  -  (i  -  x^hJct) 

+   x^Hjtt)    -   I?(Tlf,0f) 

+  El(Tl,ref'P)   "  ^^.ref'Pl.ref* 
+  Hl<Tl,ref>    "  %<Tl,ref'Pl,ref<V1,re£))] 
+   (1  -  x^)    [terms  for   component  2    similar   to   those  for   component  1] 

+  xlHl(Tl,ref'Pl,ref<Vl,ref>> 

+   <!"  xl>H2<Tl,ref'P2,ref<V2,ref>>  <4'29> 

In  this  expression,    the  first  and   second  lines   are   identical    to   those   in  Eqn. 
4.28;    the    third  line    is   the   integral    of   the   perfect    gas  molar  heat   capacity   at 
constant   pressure   between  the    temperature   of    interest  and  the   reference 
temperature   for   component  1.      The   next   line   equals  zero  because    of    the 
properties  of    the    perfect   gas.      The   fifth  line   can  be   evaluated  using  Table 
4.1.      A  similar    set   of   terms  exists  for    the   second  component.      The   final    two 
terms  on  the   right  hand   side   of  Eqn.   4.29   represent   the   pure    components   in 
their   respective   reference    states.      The   enthalpy   of    the  mixture    then  becomes 
the  following: 
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Hm(T,p(V),Xl)    =   H^T^V), xx)    -   H£(T,Xl) 


T 

.0 


+   Xl[    /      Cpl(T)    dT       +  Bl°(T1<ref)    -      ^(Tltiet,p1>tet(Vlttef)n 


Tl,ref 


+   (1  -   x^)      [terms  for   component  2    similar   to   those    for   component  1] 

+  ZlHl^l.ref'Pl.ref^ref^ 

+   (1  -   ^1)B2(T2Tef,v2fTef(\2ief))  (4.30) 

The  last    two   terms   cannot   be   determined  absolutely.      Setting   such   terms   equal 
to  zero   is  an  arbitrary   assignment   of    the   reference    state.      More   properly, 
these    terms    should  be  moved   to   the   left-hand   side    of    the   equation. 

The   entropy    for    the  mixture    is  evaluated  as  follows: 

Sm(T,p(V),Xl)   =    Sm(T,p(V),x1)   -    S°(T,V,Xl) 

+   S2(T,V,Xl)    -  XlSf(T,V)    -    (1  -  X2)S£(T,V) 

+  Xl[S?(T,V)   -    S°(Tlref,V) 

+   Sl°<Tl,ref'V>    "   Sl<Tl,ref'Vi,ref> 

(continued  on  next   page) 
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-    Sl(Tl,ref'Pl,ref<Vl,ref>>    "    S1°<T1 ,  ref » Vl ,  ref  >  > 


+  (1  -  x,)    [terms  for   component  2    similar   to   those   for   component   1] 


+  ilSl^Tl,ref'Pl,ref<Vl,ref)> 


+   (1  -   *2>S2(T2,ref'P2,ref<V2,ref>>  (4'31) 


The    path   that   is  followed   in  the    calculation  is   the    same   as  for    the   enthalpy. 
The   first,    and  fifth   lines   can  be   evaluated  from  Table  4.1   by   substituting   the 
appropriate  values  for    the  molecular   parameters   into   the   equation  of    state. 

The    second  line    can  be   evaluated  from  Eqn.   4.26.      The    third  line    is  evaluated 

o 

by    integrating  Cy/T  over   the   temperature   range   T     £    to  T.      The   fourth  line    is 

evaluated  by    integrating   (dS/dV)™  over   the  volume   range    and  using  the  Maxwell 
relation  that    OS/3V)T  =    OP/3T)y,    which   equals  R/V  for    a  perfect   gas. 
Equation  4.31    then  becomes   the  following: 

Sffi(T,p(V),Xl)    =   S^T.pfV),:^)    -    S°(T,V,x1) 

-  R(xx   n(xx)    +   (1  -  x1)£n(l  -  xx) ) 

T 
+  x^f  £v  dT  +   E  ln(V/V1   re£)) 

Tl,ref  T 

"      *<*..£■    Pl,ref    <Vl,ref>>      "   S°<Tref    Vref>l 
+   (1  -  x1)    [terms  for   component  2    similar   to   those   for   component  1] 

+  *1   h    (Tl,ref    Pl,ref   <Vl,ref>> 

+   (1  -  xx)    S2    (T2>ref,    P2jiet   (V2jref))  (4.32) 
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The   evaluation  of    any   mixture  property  proceeds   by   a   similar  method.      First 
evaluate   the   difference  between  the  real   mixture  and  the   ideal-perfect    gas 
mixture.      When  the   components  are    separated,    the  derived  nature   of    the   ideal 
mixture   prescribes   the   associated   change   in  the   property.      After    the 
separation  has  been  made,    each   component   can  be    treated  as  if   it  were  a  pure 
material. 

The   additional    degree(s)   of   freedom   with  mixtures,    the   compositions,    gives 
rise   to  a    set    of   properties   that   typically  have  a    trivial   meaning   in  pure 
materials,    the   partial    molar  properties.        Earlier   in  this   section,    the 
experiment   that  led  to  the    definition  of    the   ideal   mixture  was   discussed.      Let 
us   discuss   an  experiment  performed  under   similar  conditions  for  which  the 
outcome   is    different.      Suppose   one  had  a  mixture    of    two  materials,    A  and  B,   at 
a   particular   temperature  and  pressure.      Let  us  further    suppose    that   a    small 
amount   of  B  were  added  to  the  mixture  while   keeping  at  T  and  p   constant.      If 
one  measured  the  ensuing  change    in  volume  and  attributed  that  change   to  the 
'effective'  volume    of    the  added   component   in  the  mixture,    that   effective  molar 
volume  would  be   calculated  as  follows: 

VB  =    [V(T,p,nA,nB  +  6nB)    -  V(T,  p,  nA,  nfi)  ] /6nB  (4.33) 


where    the   bold  face  V  indicates  a   total   (rather    than  a  molar)   volume,    and  n- 
refers   to   the   total    number   of    moles  of    the  respective  components.      In  the 
limit   of  adding  an  infinitesimal    amount   of  B,    the   effective  molar  volume 
described   in  Eqn.    4.33  becomes  the  value    of    the   following  partial   derivative: 


3S 


VB  =    (6V(T,p,nA,i!B)/anB)  (4.34) 

T,  p,  n^ 


This  partial   derivative   is   the   definition  of    the    'partial   molar  volume'. 
(Partial  molar  quantities  are   indicated  by  an  overbar.)      One    could  have 
evaluated  the   partial   molar  volume   of   component  A,   V»,    in  terms  of    the   partial 
derivative   of    the  volume  with  respect   to   the   amount   of  A.      One   can  define 
similar   partial   molar  quantities  for   any   extensive   property   of    the    system.      In 
this  example,    one    can  be    comforted  by   the   explanation  that   the    partial  molar 
volumes   are     effective     volumes.      Partial   molar  quantities  often  do   not   retain 
this   simple   interpretation;    thus  one    should  keep   in  mind  the  formal  mathema- 
tical  definition  of   these   quantities. 

In   the   previous   paragraph  we  have   discussed  an  operational   definition  for    the 
partial  molar  volume.      Typically,    it   is  experimentally  determined  in  a   differ- 
ent way,    by  noting  the  variation  of    the  molar  volume  with   composition.      In   the 
following   paragraphs,    the    connection  between  the  molar  volume   and   the   partial 
molar  volume  will   be   discussed  and  an  example   of   an   important   class   of 
thermodynamic   relationships  will   be   encountered. 

Because   volume    is  an  extensive   property,    more   properly   a  first    order 
homogeneous  function  of    the   amount,    one   can  write    the   following  relation 
between  the    total    volume   and  the  molar  volume   of   a  material: 

va,?,^,^)  =  nv(T,P,  J^ir) 

=  nV(T,p,Xl)  (4.35) 
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where   n  =  n-.+  iu.      The   partial    molar  volume    of    the    two   components   can  be 
calculated  as  follows: 


1    l,p,  n2 


=  V   (T.p.iij)    +  n(  —    >         '  - 


3xl/T,p  \anl/T,p, 


p   \"»1/T,p,ll2 


=  V  +   (1 


-xl)(dx1)T) 


(where   x*    =    n^/  (n^+v^)) 


3V 
V2    (T,p,Xl)    =(^" 

2/T,p,n1 


=  V   (T.p.nj)    +   n  f  —    l       '-^ 


Sxi/T.pVani/T,?,^ 


=  V 


/av\ 


(4.36a) 


(4.36b) 


One    could  expand   this   definition  to  as  many   components  as   one  wished.      The 
resulting  partial    molar  volume  would  be 


i  i   x'P'xk^i 
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where,    8  —    is   the  Kronecker   delta    (8"=!   when   i=j    and  8,.=0  when   i^j)  ,    and     x- 
is   the    set    of   independent  mole   fractions.      Recall    that   the   number  with 
independent  mole   fractions   is   one    fewer    than  the   number   of   components    since 
the    sum  of    the  mole   fractions   is   unity. 

Figure  4.1    shows   the  molar  volume    plotted  versus  mole   fraction  for    a  binary 
mixture.      Eqns  4.36a   and  b  have   an   immediate    geometric    interpretation;    V,    and 
V2    are  respectively   the  Xi=l    and  x*=0    intercepts  of    the   tangent   to  V  at   X... 
It   immediately  follows   that 

V  =  x^   +   (1  -  x±)y2  (4.37) 

which  is   consistent  with   the   notion  that   the   partial   molar  volumes   are,    in  a 
sense,    effective  volumes  for   the  components   in  the   solution.      If   one  moves  to 
mixtures  with  more    than  two   components,    the   partial  molar  volume   becomes   the 
intersection  of   a  many   dimensional   plane    tangent   to   the  V-composition   surface 
at   the    composition  of    the  mixture  with   the  respective   pure    component  axes. 

By   evaluating  the   composition  derivative   of    the  molar  volume,    another    important 
relationship   can  be   derived. 


Wt,p     1       2       *  Wt,p  *    VsVt, 


(4.38) 


41 


LU 


O 

> 

< 


0  xB,i  1 

COMPOSITION  (mole  fraction  component  B) 


Figure  4.1.  Molar  volume  of  a  mixture  at  a  fixed  temperature  and  pressure 
showing  relationship  between  pure  component  and  mixture  molar 
volumes  and  the  partial   molar  volumes   of   each   component. 
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\ 

The   difference  V<    -   V2    can  be    evaluated  by   subtracting  Eqn.    4.36b   from  Eqn. 


4.36a;    however,    that   difference    equals  (  — 

\3X 


l7T,p 


One    concludes   then  that 


avi 


x,    "1  +    (1-x,)    2^2 

aWt,p  ^9x1 


dV, 


=  0 


(4.39) 


This  relation,    one    of    a    set   of    relations   known  as    'Gibbs-Duhem'    relations,    can 
be  justified   geometrically  by   noting    that  when   the    tangent   is  rolled   along   the 
V  -   x*    curve    that   changes    in  the   right   and  left   intercepts  are   related  by 
construction  of    similar    triangles  as  follows: 


^2      = 


SVl 

1-x-, 


For    the   purposes   of   phase    equilibrium    in  mixtures,    the  most    important    partial 
molar  quantity   is    the    chemical    potential,   which   can  be    defined   in  the 
following  ways: 


aA=    OE<S.V,{ni»/8nA)8-Vfni^  A 


OH(S.pW.{n1))/8nA)gfp  A 


OA(T,V,{n.})/anA)T)V#ni^  A 


(8«T.p(V),{ni})/a.A)T|P(ai<lA 


(4.40) 
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The   chemical    potential    for    the   Carnahan-Starl  ing-DeSantis   equation  of    state    is 
most   easily   evaluated  by   differentiating   the  Hemholtz   free   energy;    the 
resulting  expression   is   given  in  Table  4.1    for   a  binary  mixture. 

Whenever    two   or  more    phases   are   in  equilibrium,    not   only  must   the    temperature 
and  pressures   of    all   phases  be    the    same,    but   also   the   chemical    potential    of 
each   component  must   be    equal    in  all    phases.      This  result   is  a   consequence    of 
the    thermodynamic    stability    requirement    that,    at   a   fixed   temperature   and 
pressure,    the  Gibbs  free   energy   of   a    system  will    seek  a  minimum.      How   the 
phase    equilibrium   arises   can  be    seen  in  Figure  4.2,   where    the  molar  Gibbs   free 
energies  for    the    two   phase   forms   intersect   at   XqJ    however,    between  ii    and  iy, 
the  mixture   can  have  yet   a   lower  Gibbs   free   energy  by    forming  two   phases.      The 
free   energy   of    the    two   phase    system  lies  on  the    common   tangent   to   the  free 
energy   curves   for    each  phase    form.      In  an  argument    similar   to   the   one    for    the 
volumes,    the    xR  =  0    and   xB  =   1    intercepts  represent   the   respective   partial 
molar  Gibbs   free   energies   (or   chemical    potentials).      With   the   common   tangent 
construction,    one    can  readily   see    that   the   respective    chemical    potentials   of 
each   of    the   components  will   be    the    same    in  both  phases.      Were   one    dealing  with 
a    three    component    system,    there  would  be   a   common   tangent   plane    construction^ 
for    systems  with  more    than  three   components,    there  would  be    a   common  tangent 
hyperplane,    a   problem   that   is  more   easily   contemplated  mathematically   than 
visual  ized. 

In  this    section,    we  have   reviewed   several    important    ideas.      First,    we   found 
that    the    properties   of    the   perfect    gas    state   for  both  pure  materials  and  for 
mixtures  was   an  essential    device    for   connecting   one    state    of    a  material  with 
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Figure  4.2      Molar  Gibbs   free   energy   of    a  mixture   at   a   fixed   temperature   and 
pressure.      The    stable    phase    is   the   one  having   the   lower  Gibbs 
free   energy;    for   compositions  lying  between  a   common  tangent    to 
the   liquid  and  vapor   branches  a    two-phase  mixture    posses  a 
lower   free   energy   than  either  phase. 
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another.      Second,    we   found   that    the   p-V-T  equation  of    state   alone  was  not 
adequate    to   describe    the    temperature   dependence   of   a  material;    knowledge   of 
the   perfect   gas  heat   capacities  was   essential.      Finally,    we   found  that   the 
properties   of   a  mixture    could  be    connected  to   the   properties  of    its   consti- 
tuent  components   through   the    ideal-perfect    gas  mixture   and  that,    by  using  this 
device,    the   reference    states   of    the    constituent   components   of   a  mixture   define 
the  values   of    the    thermodymamic   properties  of    the  mixture    itself. 


46 


5.      CRITICAL  BEHAVIOR  AND  THE  EFFECT  OF  NEARBY  CRITICAL   POINTS 
Through  the   past  few    sections,    we  have    developed  a   fluid  model   based  upon  a 
hard  sphere   reference    fluid.      By  using  this  model,    the   thermodynamic  proper- 
ties  of   fluids,    both  pure    components  and  mixtures,    can  be    described   by  having 
only   a  modest    set  of   experimental    information.      In  this  development  we  have 
avoided  critical    points,    one   phenomenon  that    cannot   be    described  quantita- 
tively by    this  or    any   of    the  commonly  used  industrial    equations  of    state. 

First,    let   us  note    that   this  model  has  a    critical   point.      A   single    component 
fluid  is  at   its  critical    point  when  the  first   and   second  derivatives  of 
pressure  with  respect   to  volume   are  zero: 


cLe\     =   0  (&A     =   0  (5'1} 


av/T  \av2/T 


For    the  Carnahan-Starl  ing-DeSantis  equation  of    state    this  yields: 


(ta\m  o  = SI T-v4  -  bv3  -  ^£  +  iiv  -  blJ 

{dW{  v2/v       b\4L  4  "  256. 


riv-j 


+   ai2Vjt_bl_  (5  2) 


V2(V  +  b)2 


e2p\_   ft  ...  RT  Lv5    .    IbV*    .    Sb2V3   _  5b3V2    .    5b4V  _  b£. 


|=   0  = 


^%  v3 


(v  -  j; 


.^5   +  TbVl  +  5t£v 
L  2  4 


2         '  «/  v\5  2  4  16  128  512 


+   a(-6V2  -  6Vb  -  2b2)  (5#3) 

V3(V  +  b)3 


where   all   quantities  are   evaluated  at   the   critical    point. 
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Equations    (5.2)    and    (5.3)    could    in  principle   be    solved   to  yield   the    critical 
temperature   and  volume   in  terms   of    the   a   and  b   parameters*,    the    critical 
pressure  would  then  be    given  directly  by   the  equation  of    state.      The 
complexity   of    these   expressions,    however,   makes  an  analytical    solution 
impossible.      Thus  we   postulate: 


Vc  =  X1bc  (5.4) 


where  X,    is   a   constant.      Such   a  relationship  holds   for    the  van  der  Waals- 
like    equations   of    state.      Substitution  of   Eqn.    (5.4)    into  Eqn.    (5.2)    yields 
for    the  critical    temperature: 


_   ac^2  ..   ac^2  /*  o 

T«  "  m7  "  E£b7  (5-5) 


where  Xn    is  a   complex  function  of   X,    and  thus  also  a   constant.      Equation  (5.3) 

is    satisfied  using   the  above   expressions  for  V.   and  T_,    thus    confirming   the 

assumed  relationship  for  V      (Eqn.    5.4).      A  numerical    solution  of    the   above 

c 

system   of   equations  yields: 


V„  =  3.006818  b„  (5.6) 


0.227329ac 

T„  =  e  (5.7) 

Rbc 


The   critical    pressure   is   given  by 
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0.315714  RTC 
Pc  -  y (5*8> 


The   constant  0.316    in  Eqn.    (5.8)    is   the  value    of    the   compressibility   factor    at 
the   critical   point.      This   compares  to  a   critical    compressibility   of  0.375    for 
the  van  der  Waals  equation  and  1/3    for    the  Redl  ich-Kwong-Soave   equation  of 
state.      A  wide  variety   of   organic   fluids  have   a  Z_   of   approximately  0.27-0.29 
(Reid,    et  al.,    1977)  ,    thus  none    of    these   equations  of    state   can  quantitatively 
predict   near   critical  behavior.      Were  we   to  use   data  far  from   the   critical 
point    (Sengers  et  al.,    1981)    to  evaluate    the  parameters   in  the   equation,    we 
would  find  the   following  discrepencies  between  the  measured  and  predicted 
critical    properties:      first,    the  predicted  and  measured  critical    points  would 
not   coincide  within  experimental    uncertainty;    second,    the   predicted  values   of 
all    the   extensive   properties — volumes,    enthalpies,    entropies,    etc.— would 
differ  from   the  measured  values   in  fundamental  ways;   finally,    the  values   of 
the   thermodynamic  response    functions— -heat  capacities  at  both  constant   pres- 
sure  and  vol  Time,    the   isothermal    compressibility,    and   the    thermal   expansion 
coefficient — would  all   diverge  more    strongly  near   the   critical    point   than 
those  respective   properties   predicted  by   the   equation  of    state. 

One    is   tempted  to  resolve   these    differences  by   forcing  the   equation  of    state 
to  match   the    critical   behavior.      Such  coercion  does   not  represent  an 
acceptable    solution  to  what   is  a   fundamental    physical    problem.      Forcing   the 
equation  of    state   to  match   the    critical   behavior   affects   the    temperature 
dependence    of    the  molecular  parameters  a   and  b;    thus,    states   far   from   the 
critical   density — and  critical    point — but   near   the    critical    temperature  would 
'sense'    the   critical    point   in  a   physically  unrealistic  way.      The   alteration  of 
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an  equation  of    state   to  produce   both  proper   near-critical    and  far-from- 
critical   behavior   is   a  major    task  (Woolley,     1983)    (Fox,    1983)   and  will   not   be 
discussed  here.       As   long  as   one    is  not   operating   too  near   the   critical    point, 
a   somewhat  subjective   criteria   that    depends   upon   the    property  being 
considered,     'classical'    equations   of    state,     such    as   the   one   discussed  in  this 
paper,    can  describe    the   properties   of   a    fluid  quite   accurately. 

The   criteria  for   critical    points   in  binary  mixtures  are   the   following: 


D  = 


3^A 
3  V2 


T,x 


3jfA_ 
3x3V 


(Hi) 

\3x3V/ 


d±k 
3x2 


V,T 


-  0 


(5.9a) 


and 


D*  = 


3D 
3V 


3D 
3x 


T,x 


T,V 


3_2a\ 
3x3V/ 


32A 
3x2 


V,T 


=  0 


(5.9b) 


or,    by    straightforward   thermodynamic    transformation 


32G 
3x 


=    0    ®M  [  ^N)  -    0 


T.P 


3x- 


(5.10) 


T,p 


50 


Although   the   criteria  for  mixtures  appear   to  be    different   than  for   pure 
materials,    the   physical   and  thermodynamic  reasons  for    them   are    the    same.      The 
differences  are   only   the   consequences  of    the  additional    degrees   of   freedom,    the 
compositions,    available    in  mixtures  and  not   in  pure  materials.      Indeed,    the 
pure  material    criteria   expressed   in  Eqn.    5.1,    are   a    special    case    of    the 
criteria   in  Eqn.   5.9a  and  b.      The   essence    of  Eqns.    5.1,   5.9a  and  5.9b   is   the 
thermodynamic   stability   requirement   arising  from   the    second  law;    for   a    system 
to  be    stable   any   perturbation  of    the    system  at  fixed  volume  and   temperature 
must   cause    the  Helmholtz    free   energy   to  rise.      The   first   of    the   criteria   in 
Eqn.   5.1   and   that   in  Eqn.   5.9a   define    the   boundary  between  a   locally   stable 
and  unstable   thermodynamic   state.      A  critical    point    is  a    special    place    on  this 
boundary    that  requires   that  any   perturbation  will   lead   into  a    stable    state;    at 
other  points  on  the  boundary   there  will  be   perturbations   that  lead   into 
unstable    states.      The    important  fact   to   note  however   is   that   the   critical 
point  for   a  mixture   is  not   given  by    the   same   function  of   a   and  b  as   the  pure 
materials. 

Even  though   the   equation  of    state   cannot  quantitatively  predict    the   near 
critical   region,    applying   such  an  expression  to  mixtures   is   superior   to   the 
usual   mixing  rules.      The   simplest  mixture   approximation,    and  one    that    is   often 
applied,    is   that   of  an  ideal  mixture  where    the  mixture   property   (e.g.,   molar 
volume)    is   given  by    the  mole   fraction  weighted  average    of    the   pure   component 
v  al  ue  s : 

Vm  =  *AVA+  *BVB  <5'n> 
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Such  mixing  rules  can  work  well  when  both,  pure  components  are  well  below   or 
above    their   critical    points  and   if   no    strong   intermolecular  interactions   (such 
as  hydrogen  bonding)   are   present.      These   rules   can  collapse,    however,    when  one 
of   the   components  is  near   its   critical   point,    even   though   the   mixture   itself 
may   be    an  ordinary   fluid  (Morrison,    1985). 

Figure  5.1    shows  how   errors   can  arise  with   simple  mixing  rules.      In   this 
illustration,    which  plots  molar  liquid  volume  against   composition,    the 
temperature  is  near,    but  below,    the   critical    temperature   of   component  A  and 
well   below    the  critical    temperature  of   component  B   in  the  representative 
mixture  A/B.      The   curve  ab   is   the  locus   of    saturated  liquid  volumes.      The 
curve   af    is   the   isobar   at   the   saturation  pressure   of    component  A,      The    strong 
curvature    of  ab  and   af   near   the   pure  A   side    of    the  figure   are   indicative   of 
the   near   critical    state   of    component  A      The  ideal  mixture  approximation  would 
give   a    straight  line   between  a   and  b   or,   more    correctly  between  a  and   f.      In 
both   cases,    the    straight  line   deviates  from    the   corresponding  curves  by    about 
15    percent   in   the   mid   composition  range   for    this  example.      Of   course,    by   the 
nature   of    the  mixing  rule,    the   curve   and  line   must   coincide   at   pure  A  and  B.      A 
well-chosen  equation  of   state,    however,    would  follow    the   proper   curved  path, 
although  perhaps   deviating  from    the  actual   volume  as  the  mixture  approaches 
its    critical    point  at  high   concentrations   of   A. 

In   summary,    so  long  as   one    stays  away   from   the  critical   point   of   the  mixture, 
classical   equations  of   state   can  be  made   to  work  quite   adequately  for    the 
quantitative    evaluation  of    the  properties  of    the  mixture.      As  noted  before, 
one's  only  recourse    should  one    of    the    components   of    this  mixture   be   near   to   or 
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Figure  5.1.      Molar  volume   of    a  mixture   at  a  fixed  temperature   showing 
error   of   estimating  mixture   properties   by   a  linear 
combination  of    the  pure   component  values. 
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above    its  critical    point   is  the  use    of   an  equation  of    state;    all    schemes  used 
for   ordinary  liquids  break  down   completely  under   such  conditions.      As   indi- 
cated previously,    the   problem    is  more   fundamental    than  the   absence    of   liquid 
data  for   the  respective   component  above   its   critical   point;    the   presence   of  a 
nearby   critical    point   causes  properties   that   otherwise    change    slowly   and 
smoothly  with  changes   in  composition,    temperature,    or  pressure,    to   change 
rapidly.      As  a   consequence   any   linear   interpolation  scheme   ceases  to   serve  as 
a  quantitative  way   of  estimating  a   property. 
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6.      NUMERICAL   IMPLEMENTATION 

The  algorithms  which   express  the  p-V-T  and  thermodynamic  relationships 
presented  in  the   preceeding   sections  have   been   implemented   in  a    set    of 
computer    subroutines  written   in   standard  ASCII   FORTRAN  (FORTRAN  77).      These 
are    summarized   in  Table  6.1   and  fall    into  four    general    categories:       1)    subrou- 
tines  to   determine    the   a   and  b  parameters  and  the  mixing  coefficient,    f,    from 
saturation  data,    2)    routines  which    store   coefficients  for    the   a,    b  and  C° 
expressions  for  pure   components  and   calculate    these    parameters  as   a  function 
of    temperature   and   composition  for    a  mixture,    3)   property   routines  for   the 
calculation  of  enthalpy,    heat   capacity,    entropy,    specific  volume,    and    satura- 
tion pressure   and  4)    auxiliary   routines  which   are   referenced  by    the   other 
property  routines.      The    subroutines   are    interdependent  with  one    calling 
another.      This    section  details   the  methods   of    solution  employed  by    these 
routines.       The  required   inputs  and  resulting  outputs   for   each  routine    are 
summarized    in  Table   6.1    and   discussed   in  detail    in  Appendix  B. 

Determination  of  Equation  of   State   Parameters 

The   a   and  b  parameters   in  the   equation  of    state  must  be    determined  as  a 
function  of    temperature   for    each   of    the   pure  materials   of    interest.      Data   are 
required  for    this   determination.      Laboratory  measurements  or   original    litera- 
ture values   are   preferablej    tabular  values    (e.g.,     from    handbooks)    may    also   be 
used  but    are  less   desirable   because    considerable   data    smoothing  has   taken 
place   to   generate    the    tables.      The  most    commonly   available   and   generally  most 
reliable  vapor-liquid  equilibria   data   are    the    saturated  liquid  and  vapor 
specific  volumes  and   the    equilibrium  vapor    pressure   at   a    given  temperature. 
The   present   algorithm    is   developed  for    these    data   although   the  overall 
methodology  is   general   and   can  be  modified  to  accomodate    other   data. 
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Table  6.1.   Summary  of  Thermodynamic  Property  and  Associated  Subroutines 


Name De sc r ipt  ion/  Input s/  Ou t jguts 

Data  Fitting  Routines 
FITAB 

outputs:      a,b 


Other  Routines   Called 


determine    pure  component  a,b  parameters 
inputs:      T,    psatf6,  Vl>e,  Vye,    <op,  «£f  »y 


FITF  determine    interaction  parameters  from 

mixture    data 
inputs:      T,    x. ,    some    combination  of: 

psat,  e'    V£,  e'    Vv,  e'    xv,e 
output :      f  22 

Routine    to  Access    Stored  Parameters 

BCONST  access  eqn  state  parameters  from   data  base 

and   calculate  reference    states   for  H,S 
inputs:       code    nos.     for    pure    comp. ,    fQ.,  f, 
outputs:       (contained   in   common  blocks) 
note:      must   be    referenced  once    for    each 

mixture  before    calling  any   of 

following  routines. 

Property  Routines 


PL  I  MIT,    VIT 


BCONST,  BUBLT,  ENTROP 
ESPAR,  HCVCP,  PLIMIT, 
VIT,    ZXLSF 


BDESC  (via   common 
blocks),    BUBLT,    ENTROP, 
ESPAR,    HCVCP,    PLIMIT, 
VIT 


ESPAR,    PLIMIT,    VIT 


BUBLT  calculate   bubble   or   dew   point   pressure 

inputs:      T,    x   of    parent   phase 
outputs:      P,    x  of    incipient   phase,    V« ,    V 

ENTROP  compute    specific   entropy 

inputs:      T,    V,    x 
output:      S   (returned   as   function  value) 

HCVCP  compute   enthalpy  and/or  heat    capacity 

inputs:      T,    V,    x,    IQ  -    output   qualifier 
outputs:      H  and/or   C^   and/or   C      (as 
specified  by    IQ) 

Auxiliary  Routines   (transparent   to  user   but   required   ic   executable   element): 


ESPAR 


ESPAR 


BDESC 

ESPAR 

PLIMIT 

VIT 

ZXLSF 


Block  data  routine    containing  pure    component   data 


compute   a,    b,    C   ,    f^    as   functions  of   T,  x  from    stored   coefficients 


find  upper  and  lower   bounds  on   pressure    given  T 

specific  volume   for   liquid  or  vapor   as   function  of   T,  p 

or  equivalent  one-dimensional  minimization  routine  from 
math/  statistics   library    (required  only  with  FITF) 
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Pressure   and  liquid  and  vapor  volume   are  more   than  sufficient   to  determine   'a' 
and  'b'  and   thus   the    parameters   are    chosen  to  minimize    the   following   sum    of 
squares: 

ru.,x>  .  ^(^f+  »v(^-)2*  »pf-^7-f       «.i> 

where   the   e   and  c    subscripts  refer   to  experimental   data  and  calculated  value  s 
respectively   and  the  ta  terms  are  weighting  factors  which  can  be   adjusted  to 
match   the  reliability   of    the  various  measured  quantities.      The    calculated 
pressure   and  volumes  are   evaluated  using  the  equation  of    state    subject   to  the 
constraint   that   the  Gibbs  free  energy   of    the  liquid  and  vapor   phases   be    equal. 

The  function   T  is  at  a  minimum  where  the  partial  derivatives  with  respect   to 
'a'   and  'b'    are    zero: 

*£  =  o  =  J"*,' -  \.°yjl>£  +  .  /Vv,.-Vy,c\!Vv^ 


a  A        v2  J     da  ~y\       y2  / 

-C»  e  v,  e 


+  w   /psat,e Psat,c\  ?P_s_at^c  (6#2) 

P\       «2  /       3a 

N       Psat.e  ' 


*L  m  o  .  JH.*  '  Vl.c)  tlltc  +  w    /Vv.e-Vv^c^jlv^c 

3b  \     v\e       I    ab  v\      v*  )  *b 


+  u  feat.e Psat,c\   aPsat,c  (6#3) 

P\  2  /        3b 

X         Psat,e  ' 


Because    the   equation   of    state    is  fifth   order   in  volume,    Equations  6.3    and  6.4 
cannot   be    solved  explicitly.      To    set   up  an  iteration   in  a  and  b   the    calculated 

57 


quantities   resulting  from    small    changes   in  a   and  b   (6a   and  8b)   are   expanded  in 
terms  of  partial    derivatives;  for  example: 


3V/  c(a,b,T)              dVp   c(a,b,T) 
V£,c(a  +  8a'b  +  6b'T)   =  v&c(a*b'T)    +  aa 8a  +  Qh 8b      (6'4) 


Hie    combination  of   equations  6.2   and  6.3    with   equation  6.4    and    similar 
expressions    for  V„  „    and  P00+.  „   yields: 


(i)  , 


!L-(*?l*o\2   +  _^L_  /^M  +      *P       /8psattc\2l6a  + 
v2      \  da      /  y2        V  6a     J        2  \     da        /  J 

£,  e  v,  e  Fsat,  c  J 


Lvi 


i£_  dYLc  dYl.c  +  _^v_  aVv,c  avv,c  +        S 


da  db  v2  da  db 

e  v,  e 


sa 


p_  aPsat.c  aPsat.c 

da  db 

t,  e  J 


6b 


da 


Y^c  + 


^Psat.e  ~   Psat.c  )  aPsat,c 

'P\  2  /  '    da 

K  sat,  e 


(6.5) 


±L-  dV^c  6V£,c  +  _^y_  aVv,c  aVv,c 
u2  db  da  ,,2  db  da 


Lv£,e 


v,  e 


%  aPsat.c 

,2  3b 

'sat,  e  J 


8a  + 


v2    \  db   /      v2     v  db   y 

L  x.,  e  ¥v,e 


fa>p        /3P 


2        \     db 
psa t, e 


sat,  c 


6b 


(Continued  on  next   page) 
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=  J*L,'-*l,c)  m^  +         (*,,.  -  \,<)  "^  , 


\       v2  /3b  v  \        v  2  /3b 

X       V£,e  X        Vv,e 


/Psat.e        Psat.c]  aPsat,c  (6  6) 

P  V  p       2  /       9a 

Fsat,  e 


where   all   quantities   are    evaluated  at  (a,b).      The   partial   derivatives   are 
approximated  numerically,    for   example: 


*V/.c  ~  Vfc<*  +  Aa'b>   -  Vfc<*.*>>  (6  7) 

da  Aa 


3V/  c  „  V/(a,b  +  Ab)   -  V/(a,b) 

—^5  =  JL L (6.8) 

3b  Ab 


Equations   6.5    and   6.6    form    a  linear    system    in  6a   and  6b.       Given   starting 
values  of  a   and  b   the    system    can  be    solved   to    give    improved   guesses: 


,(1+1)    «   a(D    + 


6a  (6.9) 


b(i+l)    =  b(D    +  6b  (6.10) 

where    the    subscript  is   the   iteration  index. 

In  the   implementation  of    the  above  method,    the  routine  FITAB   generates 
starting   guesses   of  a  and  b  using   the    saturated  volume   data   alone.      Starting 
at  b  =   1/2  V»        the  value    of    the  a   parameter   is   expressed  in  terms   of  b  by 
equating   the   pressure    of    the  liquid  and  vapor: 
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RT(1  +  yi  +  ji2  -  Ji3)       RT(1  +  yv  +  yv2  -  yv3) 


a= (6.11) 


V  V£  +  b)        Vv(Vv  +  b) 


where 


b  .       b 

Yl   "  4V  ;         7yr  ~  4V 

I,  e  v, e 


A  Newton's  method   iteration  is   carried  out   to  find  the  value    of  b  which 
satisfies : 

Y(b)   =  0  =  G£(a,b,T,V^e)   -  Gv(a,b,T,Vv>e)  (6.12) 

If  only  specific  volume   data  were   available    this   iteration   could  be   used  to 
determine    a   and  b. 

With    the   above    starting  values  for   a   and  b   the   iteration  outlined  in  eqns. 
6.5   to  6.10  is  then  carried  out  by  FITAB.      For  each  guess  of  (a,b)   an  inner 
iterative   loop   is   used  to  calculate    the   pressure  and  liquid  and  vapor  volume. 
This   iteration  is   based  on  the   equality   of    the   liquid  and  vapor  Gibbs  free 
energy    and   is   identical    to  the   one   used   in  the  BUBLT  routine    discussed  below. 
The    iteration   is   repeated  at    (a   +  Aa,b)    and    (a,b   +  Ab)    in  order   to  evaluate 
the   partial    derivatives   (where  Aa  =   0.001    a    and  Ab  =  0.001  b).      The   system   of 
equations   6.5    and   6.6    is   then  evaluated   and    solved  for  5a  and  6b   to   generate   a 
new  (a,b)  by  eqns.  6.9,  6.10.      The  iteration  is  repeated  until  'a'  and  'b'  are 
determined   to   a    specified  accuracy. 
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The   a  and  b  determined  for    individual    data   sets  are  then  fit   to  simple 
functions  of   temperature    (such  as   eqns.    6.14,    6.15)   to   arrive   at    the   final 
expressions   required  by    the  property   routines.      Users  are  warned  that  because 
the  FITAB  routine  assumes   saturation  data   it   cannot   be   used  at   or  above    the 
critical    point.      Furthermore,    attempts  to  force    the  equation  of    state    to  fit 
liquid-vapor   equilibrium   data   near   the   critical    point   (i.e.,    at  a  reduced 
temperature  above   about  0.95)   will    cause   a   degradation  in  accuracy   for  points 
near   the  critical    temperature  but  removed  from   the   critical   region.      This 
effect   occurs  with   any   equation  of    state   that  arises  from  the   so-called     mean 
field     approximation.      This   class  includes  all    industrial    equations   of    state. 

The    subroutine  FITAB   computes  the  a  and  b  parameters  at  a  single  temperature 
from   saturation  data  at   that   temperature.      If    the  various   saturation 
quantities  are  not  all    measured  at  the   same   temperature  or  if   superheated 
vapor   data   is  to  be   used  in  computing  a  and  b,    the   subroutine  FITAB   cannot   be 
used,    rather   it   is  necessary   to   carry   out  a  non- linear  least   squares 
regression  of   the   entire   experimental    data    set.      Non-linear  regression 
routines  are   commonly  available   in  statistics   libraries  but  because   they,    as 
well   as   the   data  to  be  fit,    vary  widely   in  form  only  a   general    discussion  of 
using   such   a  routine    is   given  here.      A  regression  routine  would  find  the    set 
of   parameters  that  minimizes   the    sum   of    squares   of  residuals  between  the 
experimental    data  points  and  corresponding  calculated  values.     Typically,    an 
external      subroutine   computes  the  residuals  as   a  function  of    the   parameters  to 
be   optimized  (e.g.,   aQ,   a^,    aj ,  bQ,  b^,   b2  in  eqns.  6.14    and   6.15).     The  process 
is  iterative  and  generally  requires  reasonably   good  initial    guesses  for   the 
parameters.       For    each    set   of   parameter   estimates   generated  by    the   iteration  in 
the  regression  program  the  external   model    subroutine   is   called  to   compute  a 
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revised  set  of   residuals.      This  external    subroutine   must   therefore   update    the 
common  block  containing  the   a.    and   bi   coefficients  and  call    the  routine  BCONST 
(which   computes  reference    state  values  as  described  below)  before  computing 
the   residuals  with  BEBLT  or    similar  routines. 

Determination  of   Interaction  Parameter 

The   mixing   rule   for    the  'a'   parameter   for   a  mixture   involves  the   interaction 
parameter,    fij,    which  must  be   determined  from  experimental   data.      The  approach 
is  again  to  find  the  value    of   fj»    which  minimizes  the   sum   of    squares  of 
relative   deviations  between  measured  and   calculated  quantities: 


+  "x   <xv,e"  Vc>2  (6-13) 


Because   the  vapor   composition  must  be  between  zero  and  unity,    the  last   term   in 
eqn.    6.13    is   expressed  as  an  absolute   error.      Only   a   single  measured  quantity 
in  addition  to  T  and  x»   is  necessary   to   determine    fj^*      If   all    of    the   data 
indicated   in  eqn.    6.13    is  not   available,     the   corresponding  weighting  factors 
in   the   expression  for   T   are    set    to  zero. 

The    subroutine   FITF  will   accept   any   combination  of   Psat>   \  >   Vy  and  x^   data 
but  requires   the   a   and  b   parameters  for    the   pure    components   (i.e.,    pure 
component   data   in  the  BDESC  routine   described  below).      For   given  values  of   T 
and  %n   and   the   current   iterative  value    of   f12'     tne  BUBLT  subroutine    (discussed 
below)    is  used  to  provide    the   calculated  pressure,    volumes  and  vapor 
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composition  needed   to  compute   eqn.    6.13.      The  actual    minimization  is   carried 
out  by   the   subroutine  ZXLSF  contained  in  the   proprietary   IMSL  library   accessed 
through   the  MBS  central    computer.      It   is  expected  that  users  will   have  access 
to  a   similar  one-dimensional  minimization  routine  and  can  modify   the 
appropriate    lines   of  FITF. 

The  values  of  f,2   calculated  for  different  individual   data   sets  should  be 
constant,    or  at  most   a  function  of    temperature,    for   a   given  pair    of  pure 
components.      A  wide   variation  in  f .«    with   composition  indicates  either   poor 
mixture  data   or  an  inability   of    the   equation  of    state   to  represent    the   entire 
composition  range   of    the  mixture.      In  the   latter   case,    the   user  must    choose    a 
value   of  f|«   consistent  with  the   composition  range   of   interest.      The  values   of 
f]9  determined  near   the  pure  component  compositions  typically  will   have  grea- 
ter uncertainty   than  values   determined  in  the  middle    of    the    composition  range. 

Calculation  of  Equation  of   State  Parameters   from  Stored  Coefficients 

o 
The   coefficients  of    the   curve  fits  to  the  parameters  a,    b  and  C     for   11   pure 

refrigerants  are   stored  in  common  blocks   initialized   in  the  BLOCK  DATA  element 

BDESC      This  element  also  contains  mol  ecul  ar  weights,    critical    properties, 

reference   state   conditions  for  enthalpy  and  entropy   and  character  variable 

representations  of    the  refrigerant   names.      The  refrigerants   currently   in  the 

data  base   are   listed   in  Table  6.2.      The    data   arrays   are   dimensioned  for  20 

components   so   that  users   can  add  data   for   other  pure  materials. 

The  information  in  BDESC  for    the   particular  refrigerant   pair   of    interest   is 
accessed  by    the   subroutine  BCONST  and   stored  in   subsidiary   common  blocks  which 
are  referenced  as  needed  by   other   subroutines.      The   enthalpy  and  entropy 


63 


Table  6.2.      Refrigerants  Currently   Included   in  Data  Base 


Code  No. 

ASHRAE 

Designation 

1 

Rll 

2 

R12 

3 

P13 

4 

R13B1 

5 

R14 

6 

R22 

7 

R23 

8 

R113 

9 

R114 

10 

R142b 

11 

R152a 

Chemical  Formula 
trichlorof  luorome  thane 
di  chl  orodif  1  uor  ome  thane 
chlorotr  if  luorome  thane 
br  omo  t  r  i  f  1  uor  ome  tha  ne 
tetraf  luorome  thane 
chl  or  odif  1  uor  ome  thane 
tr if 1 uor ome thane 
1, 1 ,2-  tri  chl  oro  tr  if  1  ur  oe  thane 
1 , 2-di  chl or  ote  tr  af 1 uor  oe  thane 
1-chloro-l,  1-dif  1  uor  oe  thane 
1, 1-dif  luoroe  thane 
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routines  are   called  at   the   reference    temperatures  of   each  component   to 
establish  reference  values   of  E  and   S.       The    inputs   to  BCONST  consist    of    the 
code    numbers   for    the   pure   components   (as  listed   in  Table  6.2)    and  the   interac- 
tion coefficient  for    the  mixture    (expressed   in  terms   of    f^   and   f,    as    given  in 
eqn.    6.17).      For    a  pure   component   the   code    numbers  will  be    identical    and  f« 
and   fi    should  be    specified  as  zero.       It   is   necessary   to   call  BCONST  only   once 
(before   any    other    property    routines   are   referenced)    for    each   pair   of    pure 
components    considered. 

The   thermodynamic   property   routines   reference    the   subroutine  ESPAR  when  it   is 
necessary  to   calculate    the  a   and  b  parameters  as   a  function  of    temperature  and 

composition.      When  called  by    the   enthalpy,    entropy   or  heat   capacity   routines 

o 
the   C     for    the   pure      components  and  the   temperature   derivatives   of   a   and  b  are 

also   calculated.      This  arrangement    completely   isolates   in  ESPAR  the 

temperature   and   composition  dependence   of    the  equation  of   state   parameters 

except  for    the  Gibbs  free  energy   and  chemical    potential    statements    in  BUBLT. 

Thus,    alternative   temperature  dependencies  or  mixing  rules  can  be   accommodated 

by   changing  only    two   subroutines    (plus   of    course    the    corresponding   data   in 

BDESC)  . 

As   presented   in  Appendix  B    the  ESPAR  routine    calculates   the   a   and  b 
parameters  as: 

a   =    aQ    exp(axT  +    ajT2)  (6.14) 

b  =  b0   +  bjT  +  b2T2  (6.15) 
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The  ideal  gas  heat  capacity  is  represented  as: 


C*  =  cQ  +  cxT  +  CjT2  (6.16) 


The   interaction  parameter   is  expressed  as: 


f12   "    f0   +   flT  (6'17> 

For  mixtures  with  an  interaction   parameter   independent   of   temperature,    fi    is 
set  to  z ero. 


Values  for    the  aA,    b^   and  c^   are   tabulated  in  Appendix  A  for  11    refrigerants. 
The  value   of   the  interaction  parameter,    fio»    f°r   several  mixtures  for  which 
p-V-T-x  data  are  available  are  also  given  in  Appendix  A. 

Dew  and  Bubble  Point  Pressures 

The  hard  sphere   equation  of    state  represents  the  p-V-T  behav  ior   as  an  explicit 
expression  for   pressure   in  terms   of    specific  volume  and  temperature.      The 
application  of    the   equation  of    state   for   a  pure   component,    however,    often 
requires    saturation  pressure  as  a  function  of   temperature   alone.      For   a  mix- 
ture  the  bubble   or   dew   point   pressure   is  a   function  of    temperature  and  the 
composition  of   the  appropriate   phase  as  well.      This    calculation  is   carried  out 
by    the    subroutine  BUBLT:      given  the   temperature   of    a  mixture   (or   pure 
component)  and  the    composition  of   one   phase    this  routine    calculates   the    sat- 
uration pressure,    the   composition  of    the   other  phase,    and  the   specific  volume 
of    each    phase. 
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Vapor-liquid  equilibria  becomes  meaningless,    of   course,    above  the  critical 
temperature   and  BUBLT  returns   a  warning  message  without   carrying   out  any 
calculations  when  the   temperature  approaches  the  critical    temperature.      For 
pure    components,    a  warning   is  returned  at   temperatures   greater    than  0.99    of 
the   critical    temperature   predicted  by   eqn,    5.7.      For  mixtures   the   limit   is 
0.99   of   the   critical    temperature    of   a   hypothetical    pure  material   having   the 
same    a(T)    and   b(T)    as   the  mixture.       The   critical    temperature   of    a  mixture    is 
always  higher   than   that   of    the    corresponding  pseudo-pure  material    and   thus, 
BUBLT  will    return  a  message    of      'critical    point   exceeded'      for    temperatures 
which  may  be    significantly  below   the  mixture    critical    temperature.      This  limit 
does  not   arise   from   the  equation  of    state   itself,    but  results  from  a  failure 
of   the   solution  logic  above   the  pseudo-pure  material    critical    point.      An 
alternate  version  of  BUBLT  which  allows  computation  to  the  mixture  critical 
point  is   under   development. 

All   physically  meaningful    solutions  to  the  equation  of    state  must  have  a 
negative    slope   of  pressure  with  respect  to  volume.      Furthermore,    for    there   to 
be   both   liquid  and  vapor   solutions  for   a  pure  component  the  pressure  will    lie 
between  the   limits  where   (6p/dV)T  =  0.      These   limits  are  a  function  of 
temperature   and  are   indicated  as   Piow  and  p        in  Figure  6.1.      The    slope    of 
6p/3V  was   given  by   eqn.    5.2    and   is   repeated  here: 


(*J>)    =  RT  /y*   _  bV3   _  b¥     +  bjv  __bl\+   a(2V  +  b) 

V9V/  0/         .x4   \  4  16        256/     v2,v  +  ^ 


ylL_  b\4    V  4  16        256/      V2(V+b)2 


(6.18) 


In  the   limit   as  V    approaches  b/4   (i.e. ,   as   the  volume  approaches  the  molecular 
volume)    Op/8V)j  approaches  negative   infinity.      At  a  volume  of  V  =  3.007b  the 
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slope   is  positive  for   all    temperatures  below   the  critical    temperature  at  least 

for  values   of   a   and  b  which   are  reasonable  for  refrigerants.      (The    constant 

3.007    is   the   ratio   of    the   critical    volume   to  b(T  J.)      At   the   other   extreme,    as 

c 

V   approaches   infinity   Op/dV)™  approaches  zero  but   remains  negative.      Thus 
between  the  volumes  of  b/4   and  infinity   the   slope   dp/dV  must   twice   pass 
through  zero,    corresponding   to   Pi  ow   and  ?„„•      la   the    critical    region  of   a 
mixture   the   pressure   is  not  bounded  by   Vicm   and  p       and  the   solution  logic 
fails;    this  is  the  reason  that  BUBLT  is   limited  by   the    critical    temperature    of 
the    corresponding   pseudo-pure  material. 

The   pressure  limits  Pi^  and  p        are  found  for  a   given  a,b  and  T  by   the 
auxiliary    subroutine   PLIMIT.      Starting  at  V  =  3.007b  the  volume   is  decreased 
and   eqn.  6.18    evaluated  until    a   negative    slope   is  found.     Using   these   limits, 
the  bisection  technique    is  applied  to  volume   to  find  the  point  at  which 
Op/dV)y  =  0.      The  value   of   pj         is  then  evaluated  from   the  volume  by   eqn. 
2.6.      If   pi         is   negative,    a    small    positive  value   is   used.     To  find   the   upper 
bound  on  pressure,    the  volume    is   increased  from  3.007b  until    a  negative   slope 
is   found*,  the    bisection  method   is    then   used  once   more. 

For    a  pure   component,    the   equil  ibrium   criteria   is   the  equal  ity  of  Gibbs  free 
energy  between   the   phases: 

G£(T,V£)    =    Gy(T,Vv)  (6.19) 

The   specific  volumes  Vp    and  V      are   in  turn  functions  of  T  and  p.      Graphically, 
eqn.   6.19    is  equivalent   to   the   requirement   that   the   areas  between  a   line    of 
constant   pressure  and  the   p  vs  V  line  be    equal    as   indicated  by   the    shaded 
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areas   in  Figure  6.1.      The   general    method  of    solution  is   to   iterate    on  pressure 
starting  at   the   limits   of   p^        and  p     .      For   each   guess   of   pressure,    the 
liquid  and  vapor   volumes  are  found  as  a  function  of  T  and  p  using  the 
subroutine  VIT   (discussed  below).      The  G  for   each  phase   is   then   computed  as   a 
function  of   T  and  V   by   eqn,  4.8   (except    that   the  G      terms  will    always   cancel 
and  are    thus   omitted).      The   iteration  is   continued  by   a    combination   of    the 
secant  and  reguli-falsi  methods  until    the  computed  liquid  and  vapor  Gibbs  free 
energies  are   equal    to  within  a    convergence   tolerance. 

The   computation  of    saturated  volumes  as  a  function  of  pressure  requires  an 
iterative   calculation  which  is   contained  in  the   subroutine  VIT.     Newton's 
method  is  used  to   iterate   on  volume  until    the  pressure  computed  by   the 
equation  of    state  agrees  with   the   input   pressure.      The   previous   converged 
value   of   Vp   or  V      is  used  as  the  initial    guess.      To   speed  convergence   the 
iteration  is   carried  out   in  terms   of   ln(V)*,  for    the  vapor  volume   iteration  the 
pressure   is  transformed  into  logarithmic  coordinates  as  wel  1.      Convergence   is 
also  aided  by  constraining  guesses  of  vapor  volume   to  values   greater   than  the 
critical    volume*,   liquid  volumes  are   constrained  to  lie  between  b/4  and  V_. 

The  need  for  finding   the   upper  and  lower  bounds  on   the   pressure   iteration 
stems  from   a  potential    problem  with    solving  for  V(T,  p).     For   example,    if   a 
pressure   greater   than  p       is  input  to   the   iteration  for  vapor  volume   no 
solution  exists  on  the  vapor  branch  of   the  p-V  curve  and  the   iteration  wil  1 
either  fail    to   converge   or  will    converge  to  a   solution  on  the   liquid  branch. 
This  latter  case  would  result  in  apparently  equal    liquid  and  vapor  volumes  and 
consequently  equal  Gibbs  free  energies  for    the    two   phases.      This  would  be  a 
trivial    and  physically  meaningless   solution  to  eqn.   6.19.      This  is  a  problem 
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Figure  6.1.      p-V-T  behavior   of    a  pure  material    as  represented  by   the   equation 
of    state;   vapor-liquid  equilibria   criteria  represented  by 
equality   of    shaded  areas. 
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of    the    solution  logic  and  not   the  equation  of    state   itself.      By  bounding  the 
pressure   iteration  a   false    solution   is    avoided. 

For   a  mixture     the   criteria  for   liquid-vapor   equilibrium   (in  addition  to 
temperature  ard  pressure   equality)    are    the   equality   of    the    chemical    potential 
of   each   component    in  the   two  phases: 

|ial(T,x1,V1)    =   ua2(T,x2,V2)  (6.20) 

^lbl(T'xl'Vl)    =   »1b2(T'x2'V2)  (6.21) 

where   the   subscripts  refer   to  components  a  and  b   in  phases  1    and  2.     The 
temperature   of   the  mixture   and  the   composition  of    the   parent   phase   are   input 
quantities.      (The  phases  with   the  known  and  unknown  compositions  are  referred 
to  as    the    parent  and   incipient   phases  respectively.      For  a  bubble    point 
calculation  the   parent   phase    is  liquid*   a  parent   phase    of  vapor   gives   the  dew 
point    pressure.) 

Two  concentric   iteration  loops  are   employed   in  the  BUBLT  routine   to   solve  for 
the    saturation  pressure   and   the    composition  of    the   incipient   phase.      For   each 
guess  of   pressure   the  volume   and  chemical    potentials   of   the  parent  phase   are 
computed.      The   inner   iteration  loop  for    composition  is  entered  and   the 
incipient   phase   a,    b,    V  and  u's     are   computed  for    the   current   guess   of   x2. 
The   chemical   potentials  are   combined  with  the   current   guess   of  x2    to  yield   the 
function  used  with   secant/regul  i-fal  si    iteration: 

¥(x,)    =  — ^ x,  (6.23) 

1  Z„      +      Zr  L 
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where  zfl  =   Xj    exp    (|jial   -   |ia2)  (6.24) 

zb  =   (1  -  xx)    exp    (nbl  -   jib2)  (6.25) 

For    each    guess  of   pressure   the   composition  iteration  continues   until  ¥(x9)    is 
within  a   convergence   tolerance   of  zero.      Successive    guesses   of   pressure    in   the 
outer    iterative   loop  are  generated  by   the   secant/regul  i-fal  si  method  and  the 
function: 

<Mp)    =   1  -    (za  +    zb)  (6.26) 

This  function  goes  to  zero  when  the    chemical    potentials   are   equal    between  the 
phases. 

The   routine  BUBLT  applies   the   pure   component    iteration  for  Gibbs  free  energy 
to   the  mixture   to   calculate   the    saturation  pressure    of  a   pure    component  having 
the    same   properties  as   the   parent   phase    of    the  mixture.      This  mixture   is 
unstable  with  respect   to   two   phases  at   this   pseudo-pure    component   pressure   and 
thus    it  represents  a    lower  bound  on  pressure  for   a  bubble  point  calculation 
and  an  upper   bound  for   a    dew   point.       The   pseudo-pure    component    pressure    is 
used  as   the   starting  value    in  the  mixture    iteration.      The   second   guess   for 
pressure    (necessary   to    start   the    secant  method  iteration)    is    given  by: 


P  =    <za  +    zb)P  (for  tub^le   point) 

(6.27) 
p(2)    =   p<1)/(z     +   z   )       (for   dew   point) 


To    start   the   composition   iteration  an   initial    guess  of   equal    parent   and 
incipient   phase    compositions   is  made.       For    the    second  and   subsequent    guesses 
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of   pressure,    the  Xj    iteration  is   started  with    the   previous   converged  value    of 
Xa*      In  each   case    the    second   guess  value   for   composition  is    given  by: 


A2)    m  _ia (6.28) 

2  za  +   2b 


A  flow    chart   of  BUBLT  is   given  as  Figure  6.2.      This   routine    imposes  numerous 
checks  and  bounds  on  the    calculations.      Although   these    slow    the   execution  time 
of    the  routine    they  were  found  to  be   necessary    to   insure  reliable   convergence, 
particularly   near    the    critical   region.     BUBLT  has   been   checked   for   numerous 
refrigerant  mixtures  and  was  found  to  converge  quickly  and  reliably  for   a  wide 
range  of  temperatures.      Occasional   convergence   problems  were   encountered  at 
very   low   reduced  temperatures   (<   0.2   T  )   but   this  is   generally  well  below    the 

v 

freezing  point  where    the   equation  of    state   is  no  longer  valid  anyway. 

Enthalpy,    Heat  Capacity  and  Entropy 

The  calculation  of   enthalpy,    heat   capacity   and  entropy   are    straightforward, 
non-iterative   implementations  of   the  expressions   given  in  Section  4.      Molar 
specific  enthalpy   and  heat   capacities  at  constant  volume   and  pressure  are 
computed  by    the    subroutine  HCVCP.      Depending  on   the  value    of   an  input 
parameter   one,    two  or   all    three   of   these  quantities  are  determined  as  a 
function  of   temperature,    composition  and   specific  volume.      The  function 
subroutine   ENTROP   computes  molar  enthalpy  by   eqn.   4.49,    also  for   a  given  T,    x 
and  V.      The    calculations  for  E  and  S   employ    the  reference    state    calculated  by 
BCONST     for    each  pure   component. 
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Figure  6.2.      Flow   chart  for    the   calculation  of   dew    or  bubble   point 
pressure 
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Iteration  Technique 

Several    iterative   calculations  required  to  compute   thermodynamic  quantities 
employ  a   combination  of   the    secant  and  reguli-falsi  methods   to   converge    to  a 
solution.       The    equations   to  be    solved  are  arranged  to   give   a   function   ¥(0) 
which  is  zero  at   the    converged  value    of    the   iteration  variable,   6.      In   the 
secant   method  an   improved   guess  for  0    is  generated  by   assuming  a  linear 
function  ¥(0)   between  two   successive    guesses   and    solving   for   ¥(©)    =    0: 


6<i+»    =  0<i>    -  4<(0(i))  f.f    -6(AT"  (6.29) 

¥(0(l))    -  ¥(0(l  1}) 


where    the    superscript   is   the   iteration  index.      This   is  represented   in  Figure 
6.3   as  the  extrapolation  from  points  (2)  and  (3)   to  find  o'4  .      (If  the 
derivative   8^/30  can  be    expressed  analytically   the  bracketed  term    in  Eqn. 
6.29    is   replaced   by    1/  (9^/80)    to    generate   Newton's  method.)      The  most   recent 
values   of  0   and  ^(0)    are   retained  and   the    iteration    continues.       The   reguli- 
falsi   method  is   similar   except    that    the  new   guess  for  0    is  based  on  points 
having  a   positive   and  negative  value    of  ¥  : 


e<i+i>  =  e<p°s)  _  Y(e<Pos>)  e(pos)  -  a(neg) (6.30) 

¥(6(pos))    _  ^(neg)) 


This  method  is    shown  in  Figure  63    as   the   interpolation  between  points    (1)    and 
(3);      note    that   the   oldest    guess  for  0    is  not    necessarily  discarded. 

Given   starting  values  for  0  which  result   in  positive   and  negative  values   of 
f(0)    (such    as    the   pressure    iteration    starting  at   p,         and  p„_)    the   reguli-falsi 

J.  OW  ^]? 

method  bounds   the   solution  and  will   always    converge   for   a   continuous  function  ^(0) 
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Figure  6.3.      Graphical    representation  of    the    secant   and  reguli-falsi 
iterations   to    solve   for  ^(6)    =  0. 
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The    secant  method,    because    it    is  always  based  on  the  two  most  recent   guesses 
of  6,    often   converges  faster   but   can  also  fail    to   converge    or   even   diverge 
from    the   solution,      In  combining  these    two  methods,    a  new   guess  for  6    is 
generated  by   the    secant  method.      However,    the   information  necessary  for    the 
reguli-falsi   method  is  retained  so   that   if  &' 1+  '    given  by    the      secant   method 
lies   outside    the    bounds    set   by   reguli-falsi   (i.e.,    ©'Pos^    and   0^ne6')    a    new 
guess  for   Q^1     '    is    generated.      For    example   in  Figure  6.3,    the    secant   method 
extrapolation  from  (2)   and   (3)    to  (r   '    is    diverging   from    the    solution      in   this 
case    the   reguli-falsi    interpolation  between  (1)    and   (3)    is  used  to   compute   a 

(A  l\ 

new   guess,   &y     '. 
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7.      EXAMPLES   OF  THE  APPLICATION  OF  THE   CSD  EQUATION  OF  STATE 

In   this    section,    the    perturbed  hard-sphere   equation   of    state    developed   in   the 
preceeding   sections   is  applied   to  a   particular   class   of    fluids — the 
halogenated  hydrocarbon  refrigerants — and   their  mixtures.      First,    the 
thermodynamic   characteristics   of   11   pure  refrigerants  as   represented  by    the 
equation   of    state   are    compared  with  published  values.      For  most   of    these 
fluids,    tabular  values   taken  from  ASHRAE   (1981)    are   used,    however,    a   detailed 
discussion  for  R113   using   the    original   literature   data   is   also   included. 
Finally,    the   ability    of    the   equation  of    state    to  represent   the  behavior    of 
mixtures   is    considered  for    two   cases:      the   nonazeotropic  mixture   E13E1/R152a 
and  the    system  R22/R12  which  has   an  azeotrope. 

Pure   Refrigerants 

The    saturation  data    of   ASHRAE    (1981)*  were   used  to   generate    the   a   and  b 
parameters   of    the   equation   of    state   for  11   common  one-   and   two-carbon 
halogenated  hydrocarbon  refrigerants.      The    saturation  pressure   and   saturated 
liquid  and  vapor  volumes   over   a  reduced   temperature  range    of   approximately  0.6 
to  0.9   were    supplied   (with   equal   weighting)    to   the   programs   described   in 
Section  6    to  obtain   the   a   and   b  which  best  fit    the    data   at  a   particular 
temperature.      These    individual    values  were    then  fit    to   the   following  functions 
of    temperature: 

a  =    aQexpUjT  +    a2T2)  (7.1) 

b  =  bQ   +  bjT  +  b2T2  (7.2) 


♦The   data    of    the    'SI  Unit  Formulations'    were   used;    these    data   are    improved  and 
revised  from  earlier   ASHRAE   compilations   presented  in  English  Units. 
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These    temperature   dependencies  are   equivalent   to  those    suggested  by  DeSantis, 

et   al.    (1976)    except  for    the  addition  of    the  Tr    terms.      The   perfect    gas  heat 

o 
capacities,    C   .    were   obtained  from   literature  values   derived  from 

spectroscopic  measurements  and  were  fit   over   a  limited  temperature   range   to 

the  function: 

C*   =   cQ  +   CjT  +   CjT2  (7.3) 

The  values   of    the   a^,    hi    and   c-   for   each  refrigerant   are    given  in  Appendix  A. 

The   success  of    the   equation  of    state   in  representing  the  ASHRAE  data   is 
indicated  in  Table  7.1.      The    differences   between  the    tabulated  values  and 
those   predicted  by    the   equation  of    state    (using  the  a   and  b   given  by   Eqn.    7.1 
and  7.2)    are  listed  as  RMS  differences   over   the    temperature  range   indicated. 
The  RMS   differences  for    the   saturation  pressure  ranged  from  0.10   to  1.34 
percent  with   an  average   for   all  11   refrigerants   of  0.54   percent.      The   liquid 
volume   representations  ranged  from  0.01  percent   to  0.34   percent  with   an 
average    of  0.09    percent.      The  RMS   error  for    the  vapor  volume  ranged  from  0.14 
to  1.16  percent  with   an  average    of  0.50  percent.      These   differences  are   on  the 
same    order  as   the   accuracy   to  which   the    corresponding  experimental  quantities 
are  known   (for    example,    see   Mears  et  al.    (1955)). 

The    temperature    dependence    of    the    differences   in    saturation  pressures  and 
volumes   for  R12   and  R22   are    shown  in  figures  7.1    and  7.2.      The   smooth  curves 
are   a  result   of   comparing   the   equation  of    state  with  tabular   data    generated  by 
another   correlation  rather   than  with   experimental   measurements.      The   greater 
differences  at  higher   temperatures   are   primarily   a  result    of    the  approaching 
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Table  7-1.      Comparison  of   ASHRAE  and  Equation  of    State  Values  of    Saturation 
Pressure,    Saturated  Liquid  and  Vapor  Volume,    and  Enthalpy   of 
Vaporization  for  11  Pure  Refrigerants 


Temp. 

Range 

of 

RMS 

Difference 

(%)* 

Refrigerant 

Cor 

relation 

(K) 

Psat 

V* 

Vv 

ABva 

Rll 

240-420 

0.41 

0.04 

0.45 

0.87 

R12 

200-340 

0.26 

0.05 

0.29 

0.74 

R13 

180-272 

0.10 

0.02 

0.14 

0.17 

KL3B1 

200-300 

0.64 

0.07 

0.57 

1.46 

R14 

140-200 

0.18 

0.01 

0.18 

0.35 

R22 

220-330 

0.50 

0.06 

0.45 

1.13 

R23 

180-272 

1.34 

0.17 

1.16 

2.88 

R113 

240-43  0 

0.22 

0.03 

0.23 

0.61 

R114 

210-370 

0.72 

0.34 

0.70 

2.33 

R142b 

220-360 

0.92 

0.10 

0.83 

2.41 

R152a 

200-340 

0.62 

0.09 

0.53 

1.88 

Average 

0.54 

0.09 

0.50 

1.35 

♦RMS   Difference  = 

1 

N 

~  N    / 
Y  (ASHRAE  - 

Aj  V          EOS 

EOS  r 

100%) 

1/2 
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Figure  7.1      Percent   difference   between  ASI7P.AE  values   of    saturation  pressure 
and  liquid  and  vapor  volumes   for  EL2   and  values    correlated  by 
equation  of    state. 
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Figure  7.2      Percent   difference   between  ASHRAE  values  of    saturation  pressure 
and  liquid  and  vapor  volumes   for  R22   and  values   correlated  by 
equation  of    state. 
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critical    point;    there   is   also  an  increased  error    in  Eqn.    7.1    and  7.2    in 
representing  a   and  b  at   the   extremes   of    temperature   used   in   the    data 
correlations.      The   equation  of    state   consistently   overpredicts  the   pressure 
and  volumes  for  R22    (fig.   7.2).      Although  one  might  expect    the    differences   to 
be    centered  about   zero,    the   observed  behavior    is   the   result   of    simultaneously 
fitting   three  quantities  with   two   parameters  which   interact   in   the    equation  of 
state. 

A  more    stringent    test   for    the   equation  of    state    is   its   ability    to   predict 

quantities  not   used   in   the    correlation  of    the    parameters.      The   enthalpy    of 

o 

vaporization  is    such    a  quantity    (and  one    that    is   independent    of    C  )  .      The  RMS 

differences   between  the    predicted  and   tabulated  values    (given  in  Table  7.1) 
vary  between  0.17    and  2.88   percent  with   an  average    of  1.35   percent.      Although 
these    differences   are   larger    than   those   for    the    pressure   and  volumes   it   is 

important   to  note    that   in  many   cases   the   enthalpy   values   reported   in  ASHRAE 

o 
are   not  measured  cal  or  imetr  ically  but    are   derived  from   p-V-T  and  C 

P 

information. 

The   estimation  of    heat   capacity    involves   a    second  derivative   of    the   equation 
of    state   and   is    thus   a   particularly    stringent   test    of   physical    consistency. 
The   equation  of    state   has   been  integrated  over   the  volume    to   give    the 
Eelmholtz   free   energy  function,    from  which   all    the    other  functions,    including 
heat   capacity,    arise.      Associated  with   that   integration  is   a    temperature- 
dependent    constant  which    can  be    evaluated  by   knowning   the   perfect    gas  heat 
capacity   of    the   fluid(s)    of    interest.      As   a    test   of    the   equation  of    state   let 
us    compare  values   of   C     for    the    saturated  liquid   state    of   R152a  from    several 
sources.      The  values   calcualted  by    the   equation  of    state   use    the   perfect    gas 
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properties   presented  by    Chen  et   al.    (1975).      Also   considered  are   the  values  of 
CD  appearing  in   the   ASHRAE  (1981)    tables,    an  experimentally-based   data 
correlation  published  by   the  National   Engineering  Laboratory   (Cartwright, 
1981)    in  the  United  Kingdom,    and  a    single   experimental   value    of   C     measured  by 
Radermacher   (1983).      These   four    sets  of    information  are   compared  in  figure 
7.3.      The  four    sets  agree  well   at   the  high- temperature  end  of    the   range.       At 
lower   temperatures,    however,    the  ASHRAE  values  are  consistently  below    the 
experimental   information  of  NEL  and  Radermacher;    this   is  an  indication  of    the 
shortcomings   of   data   correlation  schemes  that  are  in  common  use.      The  most 
striking  feature    of  figure  7.3    is   that   the  values   of    C     predicted   by    the 
equation  of    state  compare  well  with  experimental    information  even  though  they 
have  been  evaluated  without  reference   to  any  previously  measured  liquid  heat 
capacities. 


Because    of    the   data    smoothing  and  correlation  necessary  to  generate   complete 
sets   of   properties   such  as   the   ASHRAE   tables,    the   use    of    original   literature 
data    to   generate   the  equation  of    state   coefficients  would  have  been 
preferable.      The   ASHRAE  tables,   while   not   always  based  on  the  most  recent 
refrigerant   data,    represent   an  extensive  compilation  in  a  uniform  format  and 
were    thus  a    convenient    source  for    the    correlation  of    coefficients.      The 
formulation  for   1,1,2-trif luorotrichloroethane    (R113)    given  in  ASHRAE   is  that 
of    Mastroianni,    et   al.    (1978).      The    consequences   of   basing    the    equation   of 
state   correlation  on  tabular  data  will  be    examined  for    this  fluid. 

Mastroianni,    et   al.    report   the    saturation  pressure    of  0.9999     pure  R113    at  26 
temperatures  with    an  estimated  accuracy   of  0.1   percent  for   pressure   and  0.02 
to  0.05   °F    (0.01    to  0.03K)    for    the    temperature   of    the   thermostated  bath*, 
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Figure  7.3      A  comparison  of    the  molar  heat   capacity   at   constant   pressure 
for    saturated  liquid  R152a    given  by    several    sources 
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saturated  liquid  density   at  9    temperatures  with   an  accuracy   of  0.1   percent; 
and   the   p-V-T  behavior   of    superheated  vapor   along  7    isochors  at  0.1    to  1.2 
times   the   critical    density    (yielding  42   data   points)   with   a  volume   uncertainty 
of  0.01   percent  and  temperature  and  pressure   uncertainties  as   stated  above. 
They   also  measured  the   critical    temperature  but   obtained  the   critical    pressure 
and  volume   by  extrapolation.      Using   these  measurements   they  fit   the    saturation 
pressure   to  a  6   parameter   expression  for  -dn  (p..+)    as  a   function  of   T  to 

Stt  t 

obtain  an  RMS  error   of  0.28   percent.      The  9   liquid  density  measurements  were 
represented  by   a  5    parameter   expression  with   an  RMS   error   of  0.12   percent. 
The   p-V-T  behavior   of    the    superheated  vapor  was  represented  with  an  RMS   error 
of  0.20  percent   using  the  Martin-Hou  equation  of    state    (involving  13 
parameters   plus   the    critical    temperature).      These    correlations  were    combined 
with   a  5    parameter   expression  for   the   perfect   gas  heat   capacity   to   generate   a 
complete    set   of    thermodynamic   tables  and  a   pressure-enthalpy   chart    (which   are 
reproduced   in  ASHRAE) . 

We  now   compare   the  literature   data  with   the  values  predicted  by   the  hard 
sphere   equation  of    state   using  a  and  b  parameters   based  on  ASHRAE  data.      The 
agreement  for    the   saturation  pressure   and   saturated  liquid  volume    (quantities 
which  were   used   in   the    correlation  of   a  and  b)    shown  in  figure  7.4    are 
excellent   for   reduced  temperatures  below  0.9    (439  K) .      The  RMS   error   over   the 
temperature  range    of    the    correlation   (240-430  K)    is  0.038   percent  for  liquid 
volumes   and  0.29  percent   for    saturation  pressures.      These    deviations  are   only 
slightly   greater    than  those    given  in  Table  7.1;    they   are   also   nearly   identical 
to  those    of   Mastroianni,    et  al .    over   the   same    temperature  range.      Thus   the  use 
of    tabular   data  for  R113   does  not  appear   to   degrade    the  fit    of    the   equation  of 
state   for   conditions  which   are  within  the  limits  of   both   the   correlation  which 
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generated  the   tabular  data   and   the   correlation  of    the   CSD  equation  of    state 
parameters. 

The   equation  of    state   parameters  were   generated  using   saturation  data   from 
240   to  430  K.      As    shown  in  figure  7.4    the   predicted   saturation  pressure  and 
liquid  volume    show    serious   deviations    (as   great  as  10   percent)    outside    these 
bounds.      This   behavior   is  an  indication  of    the   inability   of    the   CSD  equation 
of    state   to  quantitatively   reproduce    the   critical    region  rather   than  a 
consequence   of  merely  extrapolating  beyond   the  limits   of    the    correlation. 

A  comparison  of    experimental    and  predicted  pressures   as  functions   of 
temperature   and  volume   in  the    superheated  vapor   range   is   shown  in  figure  7.5. 
Apart  from   the   critical    region,    the   equation  of    state   exhibits   commendable 
accuracy  upon  extrapolation  to  higher   temperatures  and/or   away  from 
saturation.      Note    that   the  behavior    at   the   critical    temperature   is  accurately 
represented  as  long  as   the  actual    critical   region  is   avoided.      Even   though   the 
parameters  were  fit  using  only   saturation  data  below  43 CK  the   equation  of 
state   predicts   the   pressure    of   superheated  vapor   to  within  1.2   percent  at 
temperatures  as  high   as  532  K  for   reduced  volumes  between  1.43    and  10.0.      The 
hard   sphere   equation  of    state  represents   the   behavior   of   fluids   in  a 
physically  meaningful   way   and   thus,    with   care,    it   can  be    extrapolated  within 
reasonable   bounds   in  the  absence    of    data. 

The  correlation  discussed  up  to  this  point  was  based  on  the  saturation  data  of 
ASHRAE.  We  now  discuss  a  correlation  of  the  equation  of  state  based  directly 
on  the  data  of  Mastroianni,  et  al .  This  data  was  used  with  a  non-linear  least 
squares   technique   as   discussed   in  Section  6    to    simultaneously    generate    the 
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coefficients  aQ,    a1?    &2,   bQ,    bj,    b2   (the  ASHRAE-based  values  were   used  as   the 
initial   estimates   to    start    the   iteration).      Saturation   data   at  reduced 
temperatures  above  0.9   were   excluded  from   the   correlation  as  were    superheated 
vapor   data  with  reduced  volumes   below  1.5.      The    saturated  liquid  volumes  were 
given  a   greater  weight   than  the  pressure   data    to  reflect   the   greater   accuracy 
of    this  measurement  and   also  to   give  approximately   equal    total  weights  to  each 
of    the   three   data    types.      As   indicated  in  Table  7.2,    the   individual    a.    and  hi 
changed  by   as  much  as  31   percent  between  the    two   correlations.      The  values   of 
the  actual    a   and  b  parameters,    however,    changed  only   slightly   as   indicated  in 
Table  7.2   for    two  representative   temperatures.      With   the   correlation  based  on 
the  literature  data   the  fit   of    the    superheated  vapor    is   improved  at   the 
expense    of  a    slightly   degraded  fit  for   the    saturation  data.      The    total    sum   of 
squares  of    the  residuals   decreased  by   15    percent   for    the   correlation  based  on 
the  literature   data. 

This   comparison  of    the   CSD  equation  of    state  with   the   original    literature   data 
for  R113   demonstrates   the   power   of    this  expression  to  represent   the    p-V-T 
behavior   of   a  pure  fluid.      With   the   exception  of    the   critical    region,    it 
represents    saturated  liquid  volumes   to  better    than  0.1   percent  and    saturation 
pressures  and  the  pressures   of    superheated  vapor   on  the   order   of    one-half  of 
one   percent.      More    importantly    the   CSD  equation  of    state   achieves    these 

accuracies  for   both   liquid  and  vapor  with   a   single   expression  involving  only 

o 
six  adjustable   parameters   plus   the   three   coefficients  needed   to  represent   C   . 

This   is   in   contrast   to   the    three    separate   expressions   for  R113    presented  by 

Mastroianni    et   al .    involving  a    total    of  25   parameters. 
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Table  7.2      Comparison  of   Data   Correlations   for  R113   Based  on  Saturation  Data 

of   ASHRAE  and  Saturation  Plus   Superheated  Vapor  Data    of   Mastroianni 
et   al. 


Equation   of    State    Parameters 

aQ   (kJ   m3/kmol2) 
ax    (kJ   m3/kmol2  K) 
a0    (kJ   m3/kmol2   K2) 


Dq   (nr/kmol) 
b1    (n^/kmol   E) 
b«    (m3/kmol   K2 


a  @  T  =  300  K    (kJ   m3/kmol2) 

a  @  T  =  400  K    (kJ   m3/kmol2) 

b  @  T  =   300  K    (m3/kmol) 

b  @  T  =  400  K   (m3/kmol) 


Correlation  Based  on  Data   of 
ASHRAE 


»-3 


7332.59 

-2.20396  x  10^r 
-7.2656  0  x  10  7 
0.230713 


-1.87956    x  10_4 
-1.06114  x   10  7 

3545.8 
2703.4 
0.16478 
0.13855 


Mastroianni 

7489.98 
-2.34255  x  10_, 
-5.03625   x  10  7 

0.233100 


,-3 


-2.03232  x  10 
-0.82063   x   10 

3544.8 
2707.3 
0.16474 
0.13  868 


-4 
-7 


RMS   deviation  between  EOS 

Prediction  and  Data    of   Mastroianni   (%) 


1  iq  ,  sa  t 
(T  =  282 


0.038 


425   K,    5    points) 


0.057 


sat 
(T  =  238  -   436  K,    23   points)        0.38 


10.0,    19    points) 


0.85 


405     ■   532  K,    Vr  =  2.0  - 


0.51 
0.59 


Sum   of    squares   of    residuals: 


16.61 


14.19 


90 


Equilibrium  Properties   of   the   Mixture  R13B1/R152a 

As  an  example    of   a  mixture  which   does  not  have  an  azeotrope  we   examine    the 
R13B1/R152a   system.      There   is  little   published  information  about    this  mixture 
and   thus  measurements   of    the    saturation   properties  were   undertaken  at  NBS. 
The   experimental    details  and  measurements  are  presented  elsewhere    (Morrison 
and  Neal,    1985)    but   are   briefly   described  here   before   proceeding  to  a    discus- 
sion of    the   correlation  of    the   data  with   the   equation  of    state. 

Samples  for    this   study  were   prepared  by   distilling  measured  quantities   of   each 
of    the   components  from   a   gas  buret   into  a    stainless    steel    thermocompressor   at 
liquid  nitrogen  temperature.      The   amount    of   each   of    the    components  was   deter- 
mined in  two  ways:      first  by   using  the   temperature-pressure-volume  measure- 
ments from   the    gas   buret  and   the  values   of    the    second  virial    coefficient 
predicted  by    the   equation  of    state   described  earlier   in  this  paper    (Eqn. 
2.2);    second,    by  weighing   the    thermocompressor   after   each   successive  addi- 
tion of    the   components   (Morrison  and  Kincaid,    1983)  .      The   second  virial    coef- 
ficient  for  R152a  agreed  within  experimental    uncertainty  with   the  values 
listed  by  Dymond  and  Smith   (1969);    those  values  were   derived  from   the  measure- 
ments  of   Mears  et   al.    (1955).      Second  virial    coefficients  for  E13B1  were   not 
available   in  the  literature.      The   two  amount   determinations   typically   agreed 
with  one   another   to  0.1%.      The  measurements   included  runs  on  each   of    the   pure 
refrigerant  materials  and  on  mixtures   that  were   roughly  25,   33,   50,    67,    and  75 
mole   percent  R152a.      The  mixtures  were  moved  from   the    thermocompressor   into 
the    sample   cell,    constructed  from   a   drawn  sapphire   tube   having  a  volume   of 
approximately  7.0   cm     (Davis,    1983).      The  volume  accessible    to   the    sample 
could  be    changed  by   raising  and  lowering  the  mercury    level    in  the    sapphire 
tube. 
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The  cell  was  kept   in  a  water  bath  whose   temperature  was  measured  with  a  quartz 
crystal   thermometer   calibrated  to  an  accuracy   of  0.001   K  with  an  NBS- 
calibrated  25   ohm  platinum   resistance   thermometer*,  the  temperature  was 
controlled  to  within  0.0003   E  of  a   constant  value.      The   total  volume   acces- 
sible  to  the   sample  and  the  volumes  of   the   individual   phases  were  determine  by 
measuring  the  distances  between  the   top  of   the    cell   and  the  liquid-vapor 
meniscus   or   the  liquid  mercury  meniscus.      The  volume  of   the  cell  was  cal  i- 
brated  to  within  0.001    cm     with   triply   distilled  mercury.      Pressures  were 
measured  with   a  differential    gauge    calibrated  to  0.2   kPa  with  a  dead  weight 
gauge. 

Measurements  of    the  liquid  and  vapor  volumes  and  the  pressure  were  made   at 
five   equally   spaced   temperatures  between  15°C  and  55°C.      Several    sets  of 
measurements  were  made  by  progressively  enlarging  the  volume  accessible   to  the 
sample.      The  results  most   immediately  determined  from   these   data  are   the 
liquid  molar  volume  and  the  pressure   on  the   bubble  line.      Although  no   samples 
of    either  phase  were   taken  during  the   experiments,    there  is   sufficient   infor- 
mation in  these   data   to  locate   the    dew   point   curve    (Knobler   and   Scott,     1982). 

The  bubble   point   pressures   are    shown  as   open  circles   in  figure  7.6.      The 
filled  circles  represent   the   composition  of   the   dew   phase   in  equilibrium  with 
the  50%   mixture   calculated  from    the  data.      The  molar  volumes  of   the   saturated 
liquid   are    shown   in  figure   7.7. 

The   equation  of    state  was   generated  by    first   fitting  the   properties  of   the 
pure  materials   in  the  fashion  refered  to  in  Section  6.      A  value    of  0.089    for 
the  mixing  parameter  was  found  to  optimize  the  pressure  correlation  along 
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Figure  7.6      Saturation  pressure   for    the  mixture  R13B1/R152a.      The   open 
circles,  0,    are   experimental   bubble   points  and   the    closed 
circles,  #,    dewpoints.      The   solid  and   dashed  curves  are    the 
predicted  bubble   and  dew   curves  respectively  beginning  at    the 
lowest   pressure   and   temperature,    15.902°C  and   incresing   in 
temperature   to  25.323°C,    35.323°C,    45.134°C,    and  55.033°C. 
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Figure  7.7      Saturated  liquid  volumes  for    the  mixture  R13B1/R152a.      The   open 
circles  are    the  measured  values  and   the  lines,    the  values 
predicted  by   the   equation  of    state.      The   temperatures   are 
identical    to   those    of   Figure  7.6. 
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the  babble   curve.      We   elected  to  use    only   the   pressure   data   and  not    include 
the  volume   data   to  find   the  mixing   parameter   because    of    the  higher   precision 
of    the  pressure  measurement. 

The  pressures   and  liquid  volumes   given  by    the  equation  of    state   are    shown  in 
figures  7.6   and  7.7    as    solid  and  dotted  curves   for   the   bubble  and   dew  lines 
respectively.      The  pressures  are   correlated  to  +/-  5   kPa;    the  predicted 
volumes  for    the  mixture   are   consistently   greater    than  the  measured  volumes. 
One    should  note   the  dramatic  narrowing  of    the   two  phase   region  on  the  R13B1 
side   of    the   phase   diagram,    especially   at   the  highest   temperature.      This 
behavior   does  not   indicate   the   onset   of    an  azeotrope,    rather   it    is  an 
indication  of    the  nearness  of    the  R13B1   critical    point   at  67°C. 

Azeotropes   and   the  Equilibrium  Properties  of  R12/R22 

In  the  above   discussion,    we   examined  the    simplest   kind  of   behavior   one    can 
expect  from  a  mixture    of   two  liquids,    complete  miscibility  and   the  monotone 
variation  of   properties  from   one    component   to  the   other.      We   now   consider   an 
example   of    the   next  most    complex  behavior   case,    azeotropy.      Our  motivation  is 
twofold,    first,    to  discuss   the   situations  when  one    can  expect   an  azeotrope, 
and   second,    to  examine   a  refrigerant  mixture  for  which   there   is    calorimetric 
data. 

Let  us   consider   the   situation  when  one  would  expect   an  azeotrope.      An 
azeotrope   is   inevitable  when  the    saturation  curves   of    the    two   components   cross 
the   p-T  projection.      Such   a  point   of    apparent   intersection  is   called  a 
'Bancroft   point'    (Rowlinson  and  Swinton,    1982b).      The  mixture   R12/R152a,    also 
known  as  RSOO,    shows  just    such   a  behavior.      In  mixtures  where   there   is  a 
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Bancroft  point,    the  azeotrope    can  exist   over   the   entire   composition  range   from 
one   pure    component   to   the    other.       In  R12/R152a    the   locus   of    azeotropes   is 
intercepted  by    the    solid  phase    and   the   critical    locus    (Pennington,    1952).      A 
Bancroft   point   is  not  required  for   an  azeotrope   to  exist.      Since   any  kind  of 
non-ideal   behavior   causes   a   deviation  from   the  Raoult   law   prediction,    the 
closer    the  vapor   pressures   of    the    two   components,    the  more  likely    there  will 
be   an  extremum   in  the  T-x  or   p-x  lines  and,    hence,    an  azeotrope.      Closeness   of 
vapor   pressures   is   often  associated  with   closeness   of    critical    temperatures, 
since,    with   a  few   notable   exceptions,    the   critical    pressure   of   many  materials 
(including  halogenated  hydrocarbons)    are   nearly   the    same,    about  4   MPa,    and 
saturation  lines  are   roughly  parallel    in  a   p-T  projection.      One    can  thus 
conclude    that    two  materials  with  nearly    the    same    critical    temperature   are 
likely   to    show    azeotropic   behavior.      A  sharp  demarcation  cannot  be    drawn, 
however.      A  certain  critical    temperature   difference   needn't   guarantee    the 
presence    or   absence    of    an  azeotrope;    the   situation  depends   upon  the  detailed 
molecular   character   of    the    components   in  the  mixture. 

The  mixture  R12/R22,    also  known  as  RS01,    has   an  azeotrope    that   emerges  from 
the   pure   R22   axis   of    the   phase    diagram  at   approximately  320  E  and  moves   into 
the  mixture   region  as   the   temperature    is  lowered  until,    at  253   K,    the 
azeotropic    composition  is  approximately   10  mol    percent  R12   (Spauschus,    1962). 
The  variation  of    the  azeotropic   composition  with   temperature   is  not   unusual. 
The    amount   by  which   it  varies  with   the   temperature  will   however,    be    different 
from   one   azeotropic  mixture    system   to  another.      The  azeotrope    in  R12/R22    is   a 
positive   azeotrope    (Rowlinson  and  Swinton,    1983b);    that   is,    the   boiling  curve 
will   have   a  maximum  pressure  when  measured  at   constant   temperature   or, 
conversely,    a  minimum    temperature  when  measured  at   constant   pressure.      This   is 
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the   kind  of   azeotrope    one  would  typically   expect   from  refrigerant  mixtures. 
The   opposite   kind  of   behavior   is   usually  encountered  when  the    two   components 
have   a    strong  attractive   interaction,    such   as   the  hydrogen  bonding  encountered 
in   chloroform/acetone    (Karr  et   al. ,    1951). 

The  mixture  R12/E22    shows   azeotropic  behavior    over   a  range    of    temperatures  and 
compositions  as    shown  in  fig.   7.8.      The   azeotropic  point   is   typically  found  by 
locating  an  extremism    in  the  bubble   point   curve;    at   such   a   extremum,    the   second 
law   of    thermodynamics  requires   the  liquid  and  vapor   phases  to  have    the    same 
composition  and  the   dew   and  bubble  lines   to  be    tangent    (Bett   et   al.,    1975). 
Because   of   the   flatness   of   these  lines   around  an  azeotrope,    the   uncertainty   in 
the   composition  is   typically   large.      The  data   of   Eiseman   (1957)  ,    for   example 
determine    the   azeotrope   only   to  within  5   mole   percent  as    shown  by    the      error 
bar   in  figure  7.8.      The  azeotropic   composition  can  be    determined  more 
accurately  by  a   distillation  process    (Pennington,    1952). 

The   equation  of    state   parameters  for   pure  R12    and  R22  were   determined  by   using 
data  from  the   ASHRAE  Tables    (1981)    as   described  above.      The  mixing  parameter 
was   determined  by   fitting   the  bubble   point  predicted  by    the   equation  of    state 
to   the   constant   pressure   boiling   temperatures  measured  by  Eiseman   (1957)    over 
the  full   range   of    compositions.      The  values  of    the  mixing  parameter  were 
slightly   composition  dependent.      Equally  weighting  all  values   over    the   entire 
composition  range  yielded  an  average    of  0.041  with  a    standard  deviation  of 
0.009. 

The  first   test   of    the   equation  of    state   is   the   comparison  of    its  prediction  of 
the   temperature   dependence    of    the   azeotropic   composition.      The    solid  curve   in 
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Figure  7.8  The  variations  of  the  azeotropic  composition  of  the  R12/R22 
mixture  as  reported  by  several  authors  and  predicted  by  the 
equation  of    state. 
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figure  7.8    shows   that   this  prediction  falls  within  the   experimental 
measurements   of    the   azeotrope.      As  we  have   noted  previously,    there   is 
considerable   uncertainty    in  those   points.      The   dashed  line    in  the   figure   is 
the  locus   suggested  by   Spaschus   (1962)    from  measurements  he  made   on  a    grid  of 
composition- temperature   conditions.      Our   prediction,    which   arises  from  data   at 
a    single   pressure   nearly   coincides  with  Spauschus's  experimental   locus. 

The   equation  of    state   contains   enough   information  to  evaluate   enthalpies  as 
long  as   the    temperature   is  fixed.      (For   differences   in  enthalpy  due   to 
temperature   changes,    one   also  needs   the   perfect   gas  heat   capacities.)      Neil  son 
and  White    (1959)    have  measured   the   enthalpy   change   associated  with   the 
complete   evaporation  of   mixtures  of   R12/R22   over   the   entire   composition  range 
at  222.0  K.      Their   data,    plotted   in  figure  7.9    are   integral  quantities 
because    the   composition  of    the  liquid  changed  as   the   evaporation  proceeded. 
The   prediction  of  a    closely  related  quantity,    the   enthalpy   change   associated 
with   the   complete  vaporization  of    the  liquid  mixture  at  fixed  composition  to 
its  vapor   at   the    same  fixed  composition,    is    shown  by    the    solid  curve   in  figure 
7.9.      The  agreement  between  these    two   closely   related  quantities   is  on  the 
order   of  +/-  0.5%.      The    equation  of    state   also  predicts   nearly  quantitatively 
the   curvature   associated  with   the   composition.      By  using  these    evaporation 
data,    Neilson  and  White  were   also  able   to   evaluate   a  quantity   closely  related 
to  the   enthalpy   of  mixing  along   the   saturation  line,    He.      The   comparison  of 
their   calculated  values  and   the   predicted  values  are    shown  in  figure  7.10.      The 
equation  of    state   overpredicts   this  quantity.      One    should  note  however   the 
small  magnitude    of   this  quantity.      Furthermore,    the   experimental   and  equation 
of    state  quantities   are   not   exactly   the    same:      the   experimental   value    should 
underpredict    slightly.      This   comparison  with   calorimetric   data   is  encouraging. 
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Figure  7.9      The  heat   of    complete  vaporization  at  222. 00K  for    the  mixture 

R12/R22    as  measured  by  Neilson  and  White    (1959)    and  predicted  by 
the   equation  of    state. 
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Figure  7.10     The  heat   of   mixing  for    the  R12/E22  mixtures   as  reported  hy 

Neilson  and  White    (1959)    and  predicted  by   the   equation  of    state, 
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That    such   intricate   details  about   this  mixture  are  accessible   through   so 
little   information  is  a   demonstration  of    the   power  and  versatility   of   the 
Carnahan-Starl  ing-DeSantis   equation  of    state. 
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8.      CONCLUDING   REB1ARKS 

This  Note   presents  the  development   of    an  equation  of    state   based  on  a  hard- 
sphere  reference   fluid*,  the   expression  used   is   a  modification   of    the    Carnahan- 
Starl  ing-DeSanti  s    (1969)    hard-sphere   fluid   proposed  by   DeSantis,     et   al    (1976). 
The  resulting  expression  is   simple,    based  on  a    good  physical    model   and  has   a 
physical    interpretation  of    its  parameters.      The  complete   set  of   thermodynamic 
functions  derived  from   this  equation  of   state  represent  with   thermodynamic 
consistency    the  behavior    of   both  liquid  and  vapor   phases  of   pure   fluids   and 
binary  mixtures.      The    CSD  equation  of    state  has   been  demonstrated  to  represent 
the   p-V-T  properties  of    the  11  pure  refrigerants   investigated  very   well,    often 
within  the  accuracy  to  which  the  experimental  quantities   are   known.      It   is 
likely    to  do   equally  well    for   other    similar  materials   (e.g.,    simple  hydrocar- 
bons and  hydrocarbon  derivatives)   and   simple    inorganic   molecules. 

The   representation  of   mixture  properties  with  a  single   interaction  parameter 
(in  addition  to  pure   component   information)   was   demonstrated,    detailed  aspects 
of   mixture  behavior   were   accurately  predicted  even  though   the  correlation  of 
the   interaction  parameter  was  based  on  limited  data.      There   is,   however,    a 
wide   range    of  behaviors  possible   in  mixtures  and  this   is  an  area   that  requires 
continued  study.      The  present  mixing  rules   need  to  be   applied  to  many  more 
fluid  pairs.      In  addition,    the  modification  of    the  Carnahan-Starl  ing  expres- 
sion to  account  for   a  mixture   of  hard   spheres   of   differing  diameters   (dis- 
cussed  in   Section  3)    should   be    investigated. 

The   CSD  equation  of   state   concisely  represents   properties  with   thermodynamic 
consistency   and  permits  a   rational    estimation  or   extrapolation  of  quantities 
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for   which  no  measurements  are   available.      It   is  not,    however,    intended  to 
replace  experimental   data   or  preclude    the    need  for   experimental    data    covering 
the    entire   range    of    interest. 

The  representation  of   both  liquid  and  vapor   phases  with   a   single,    physically 
well-founded   equation  of    state   is   the    desirable    state    of    affairs.       Therefore, 
this,    or    similar,    equations  of    state    should  be    considered  as   the  basis  for    the 
preparation  of   property   tabulations  for   inclusion  in  handbooks,    etc.      Such  a 
preparation  must   be   based  on  a   correlation  of   carefully  evaluated  experimental 
data   (either   original    or  literature  values).       Most    of    the   numerical 
coefficients  presented   in  this  work  are  based  on  a   correlation  of   tabulated 
data  because   of   project   time  limitations  and   in   some    cases   the   unavailability 
of    the   original    data.      For    this  reason,     calculated  values  based  on  these 
coefficients   are   not  meant   to  replace    compilations    such   as    ASHRAE   (1981). 

At   present,    the  primary  use    of    the  equation  of    state   is  in  studies  where   the 
properties   of  a  number   of   fluids   or  mixtures   are  required.      A  single    set    of 
computerized  property   routines  would  be   usable  for  different  fluids   simply  by 
changing  the  numerical   coefficients.      Because    of    the    small   number   of    coeffi- 
cients required,     it  would  be   of   particular  use  with  fluids  for  which  limited 
data    is    available. 
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Table  A.2.      Coefficients  for   Curve-fit   of    Pur e- Component  Perfect-gas  Heat 
Ca  pa  c  i  ty 


Refrigerant 

c0 

cl 

104c2 

Temp.    Range    of 
Data   Fit    (K) 

Source 

Rll 

22.0418 

0.260895 

-2.45319 

200-400 

JANAF 

R12 

17.5387 

0.248546 

-2.16271 

200-400 

JANAF 

R13 

13.93  00 

0.232181 

-1.82929 

200-400 

JANAF 

R13B1 

19.9537 

0.216394 

-1.70241 

200-400 

JANAF 

R14 

11.0629 

0.209740 

-1.40992 

200-400 

JANAF 

R22 

17.0547 

0.161633 

-0.91256 

200-400 

JANAF 

R23 

20.4760 

0.106183 

-0.12189 

200-400 

JANAF 

R113 

76.2637 

0.119641 

0.71879 

240-420 

ASHRAE 

R114 

20.7005 

0.464035 

-4.175  89 

240-470 

ASHRAE 

R142b 

23.7611 

0.231706 

-1.06534 

250-500 

Mears,    et    al 

R152a 

22.2804 

0.154009 

-0.03067 

200-400 

Chen,    et   al. 

NOTE:      Cp(kJ/k  mol   K)    =   cQ  +   CjT  +    ^T2 

where    temperature  T  is  expressed   in  Kelvin 
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APPENDIX  B :      FORTRAN   SUBROUTINES 

This   appendix  contains  listings   of    the  FORTRAN   subroutines   developed  to 
implement    the   equation  of    state   for   binary  refrigerant  mixtures  as  well    as  a 
sample   run  which   uses   the    subroutines   to   compute  a    table   of   properties  for    an 
example  mixture. 

Notes   on  the   use    of    these    routines: 

1)  Language.      All    of    the  routines   are  written  in  ASCII    standard  FORTRAN 
(FORTRAN  77).      Generic   names   are   used  for    all    calls   to   intrinsic 
functions    (e.g.,   EXP  and  LOG) . 

2)  Argument   lists.      Inputs  and  Outputs  for    each   of    the    subroutines   are 
described   in   comment    statements  at    the   beginning   of   each  routine. 

3)  Units.      The    subroutines   as   presented  here   require   all    inputs  and  compute 
all   outputs   in  SI  units  with   the   following  multipliers: 

composition:  mole   fraction 

temperature:  K 

volume:  m°/kg  mol    (equivalent   to  £/g  mol) 

pressure:  kPa 

enthalpy:  kJ/kg  mol 

entropy,    heat  capacity    and   gas   constant:      kJ/kg  mol    K 

The  routines  will  work  with  any   other    set    of    consistent   units    if    the 
values   of    the   gas   constant   and  equation  of    state    coefficients   are 
converted   in  the  BLOCK  DATA  routine.      The    user   is  reminded,    however,    that 
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English   units    (ftJ,    psia,    Btu,     lb,    etc.)    are   not   consistent   and  would 
require   insertion  of    conversion  factors   into  each  of    the   property 
subroutines. 

4)  Convergence    tolerances  and  machine   precision.      The   convergence    tolerance 
for    the   iteration  loops   in  the    subroutines  BUBLT  and  VIT  are    set   by   the 
value    of   TOLR   in  BLOCK  DATA.      The  BUBLT  routine    contains  nested  iteration 
loops.     As   the    calculation  proceeds  from   the   inner   to  outermost  loop  the 
convergence    tolerance   must  be   progressively  relaxed/   this  is  accomplished 
by   scaling  the   convergence   tolerance   of   each  loop  to  a  multiple    of  TOLR. 
A  value    of   10        is  presently  used  for  TOLR;   this  value  worked  well   on  a 

Sperry-Univac  1100/82  mainframe   computer  with  a   single  precision  word  of 

—3  fi 
36  bits   (approximately  8   decimal    digits  of    accuracy)    and  a   range   of    10 

to   10     .      This  yielded  a  final    precision  of  about   one   part   in  10      for    the 

calculation  of    saturation  pressure.      With   a  machine    of   lesser   precision, 

a  larger  value   of  TOLR  may  be  required.      Greater   precision  can  be 

obtained  by   reducing  the  value    of   TOLR  but   this  may  require  converting 

the   entire    set   of   subroutines  to   double   precision.      (The   programs   are 

presented  in   single   precision  except   for    subroutine  VIT.) 

5)  Warning  messages.      The  logical   unit  to  which  any  warning  messages   are 
written  is   currently   set   to  6    in  the  BLOCK  DATA  routine;   this  may  be 
changed  to   suit   different    systems. 
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1  SUBROUTINE   FITAB    (T.PE, VLE,  WE.WP.WL.WV.A.B.PC.VLC, WC) 

2  C 

3  C  THIS  SUBROUTINE  CALCULATES  THE  *A*  AND  'B*  PARAMETERS  OF  THE 

4  C  EQUATION  OF  STATE  WHICH  BEST  FITS  A  SET  OF  EXPERIMENTAL  SATURATION 

5  C  DATA  FOR  A  PURE  COMPONENT  AT  A  GIVEN  TEMPERATURE.   THE  REQUIRED 

6  C  DATA  ARE  SATURATION  PRESSURE  AND  LIQUID  AND  VAPOR  MOLAR 

7  C  VOLUMES  AT  TEMPERATURE  T.   WEIGHTING  FACTORS  ARE  ASSOCIATED  WITH 

8  C  EACH  EXPERIMENTAL  QUANTITY  AND  MAY  BE  ADJUSTED  TO  ACCOUNT  FOR  THE 

9  C  RELATIVE  RELIABILITY  OF  THE  VARIOUS  EXPERIMENTAL  QUANTITIES.   THE 

10  C  INDIVIDUAL  A  AND  B  VALUES  GIVEN  BY  THIS  ROUTINE  WOULD  THEN  BE  FIT 

11  C  BY  A  SEPARATE  ROUTINE  TO  AN  APPROPRIATE  FUNCTION  OF  TEMPERATURE 

12  C  FOR  USE  WITH  THE  EQUATION  OF  STATE. 

13  C 

14  C  INPUTS: 

15  C  T  -  TEMPERATURE  (K) 

16  C  PE  -  EXPERIMENTAL  SATURATION  PRESSURE  (KPA) 

17  C  VLE  -  EXPERIMENTAL  LIQUID  PHASE  MOLAR  VOLUME  (M**3/KMOL) 

18  C  WE  -  EXPERIMENTAL  VAPOR  VOLUME  (M**3/KMOL) 

19  C  WP  -  WEIGHTING  FACTOR  FOR  PRESSURE  DATA 

20  C  WL  -  WEIGHTING  FACTOR  FOR  LIQUID  VOLUME  DATA 

21  C  WV  -  WEIGHTING  FACTOR  FOR  VAPOR  VOLUME  DATA 

22  C 

23  C  OUTPUTS: 

24  C  A  -  EQUATION  OF  STATE  PARAMETER  ASSOCIATED  WITH  INTERMOLECULAR 

25  C  ATTRACTION  (KJ  M**3/KMOL**2) 

26  C  B  -  EQUATION  OF  STATE  PARAMETER  ASSOCIATED  WITH  THE  MOLECULAR 

27  C  VOLUME  (M**3/KMOL) 

28  C  PC  -  SATURATION  PRESSURE  CALCULATED  BY  EQUATION  OF  STATE  USING 

29  C  ABOVE  VALUES  OF  A  AND  B  (KPA) 

30  C  VLC  -  CALCULATED  LIQUID  VOLUME  (M**3/KMOL) 

31  C  WC  -  CALCULATED  VAPOR  VOLUME  (M**3/KMOL) 

32  C 

33  C  OTHER  SUBROUTINES  REFERENCED: 

34  C  PLIMIT  -  DETERMINES  INITIAL  BOUNDS  ON  PRESSURE  ITERATION 

35  C  VIT  -  ITERATION  FOR  MOLAR  VOLUME  GIVEN  T,  P 

36  C 

37  C  NOTE:   IT  IS  NOT  NECESSARY  TO  REFERENCE  BCONST  WHEN  USING  THIS 

38  C  ROUTINE.   ALL  NECESSARY  COMMON  BLOCKS  ARE  INITIALIZED  HERE. 

39  C  THE  GAS  CONSTANT  IS  SET  (IN  SI  UNITS)  IN  THIS  ROUTINE. 

40  C 

41  IMPLICIT  REAL  (A-H.O-Z) 

42  COMMON  /TOL/  TOLR, ITMAX, LUP 

43  COMMON  /RDATA4/  R 

44  DIMENSION  AA(2) ,BB(2) ,VL(2,2) , W(2,2) ,P(2,2) ,PG(3) 

45  1   ,PL(3),FP(2) 

46  LOGICAL  LPPOS, LPNEG, LV1CON, LV2CON 

47  C 

48  C  STATEMENT  FUNCTIONS  FOR  GIBBS  FREE  ENERGY  AND  THE  DERIVATIVE 

49  C  OF  GIBBS  FREE  ENERGY  WITH  RESPECT  TO  THE  B  PARAMETER 

50  C 

51  G(T,V,A,B)=R*T*(-LOG(V)+0.25*B/(V-0.25*B)**3 

52  1   »((8.0*V-2.25*B)*V+0.1875*B»B))+A/B*LOG(V/(V+B))-A/(V+B) 

53  DGDB(T,V,A>B)=R*T/(V-0.25*B)»*3*((-2.25*V+0.375*B)+ 

54  1   ((2.0«V-0.5625*.B)*V+0.046875*B*B)*(V+0.5*B)/(V-0.25*B))+ 

55  1   A*((1.0/(V+B)-ALOG(V/(V+B))/B)/B+1 .0/(V+B)**2) 

56  R=8.314 

57  TOLR=1.0E-7 

58  ITMAX=20 

59  LUP=6 

60  PE2=PE*PE 

61  VLE2=VLE«VLE 

62  WE2=WE*WE 

63  RT=R*T 

64  C 

65  C  USE  EXPERIMENTAL  LIQUID  MOLAR  VOLUME  AS  INITIAL  GUESS  FOR  B  AND 

66  C  ENTER  ITERATION  TO  FIND  A,  B  WHICH  EXACTLY  SATISFY  VOLUME  DATA 

67  C 

68  B=VLE 

69  DO  200    IT=1 .100 

70  YL=0.25*B/VLE 
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71  YV=0.25*B/VVE 

72  A=RT*((1 .0+(1 .0+(1 .0-YL)*YL)*YL)/(VLE*(1 .0-YL)**3) 

73  1   -(1  .0+(1  .0+(1  .0-YV)*YV)*YV)/(WE*(1  .0-YV)**3))/ 

74  1   (1  .0/(VLE*(VLE+B))-1  .0/(WE*(WE+B))) 

75  FG=G(T,VLE,A.B)-G(T,WE,A,B) 

76  IF  (ABS(FG).LT.TOLR)  GOTO  240 

77  B2=B-FG/(DGDB(T,VLE,A,B)-DGDB(T,WE,A,B)) 

78  IF  (B2.GT.4.0*VLE)  B2=(B+4.0*VLE)/2.0 

79  IF  (B2.LT.0.0)  B2=0.5*B 

80  B=B2 

81  200  CONTINUE 

82  240  AA(2)=A 

83  BB(2)=B 

84  DO  260  1=2,1,-1 

85  DO  260  J=3-I,2 

86  VL(I,J)=VLE 

87  260  VV(I,J)=VVE 

88  C 

89  C  ENTER  ITERATION  TO  FIND  A,  B  WHICH  BEST  FIT  P  AND  V  DATA 

90  C 

91  DO  400  IT=1 .ITMAX 

92  C 

93  C  FIND  VOLUMES  AND  PRESSURE  PREDICTED  BY  EQUATION  OF  STATE  FOR  CURRENT 

94  C  GUESS  OF  A,  B  AND  FOR  SLIGHTLY  DIFFERENT  A,  B  (FOR  PURPOSES  OF 

95  C  NUMERICAL  DIFFERENTIATION) 

96  C 

97  AA(1)=1 .001*AA(2) 

98  BB(1)=1 .001»BB(2) 

99  DO  360  1=2,1 ,-1 

100  DO  360  J=3-I,2 

101  C 

102  C  CALL  SUBROUTINE  TO  DETERMINE  THE  UPPER  AND  LOWER  BOUNDS 

103  C  ON  PRESSURE  FOR  WHICH  THERE  ARE  BOTH  LIQUID  AND  VAPOR 

104  C  SOLUTIONS  OF  THE  EQUATION  OF  STATE 

105  C 

106  CALL  PLIMIT  (T.AA( I ) ,BB(J) .VLOW.VUP, PLOW, PUP) 

107  C 

108  C  SET  INITIAL  GUESSES  FOR  PRESSURE  NEAR  THE  UPPER  AND 

109  C  LOWER  BOUNDS.   IF  THE  LOWER  BOUND  FOR  PRESSURE  IS  NEGATIVE 

110  C  RESET  IT  TO  A  SMALL  POSITIVE  VALUE. 

111  C 

112  IF  (PLOW. LE. 0.0)  THEN 

113  VLOW=0.8*BB(J) 

114  TC=AA(I)/(BB(J)*4.398909*R) 

115  PC=0.1049995*R*TC/BB(J) 

116  PL0W=1 .0E-12*PC 

117  PG(1)=PL0W 

118  ELSE 

119  PG(1)=PLOW+0. 0001* (PUP-PLOW) 

120  END  IF 

121  PG ( 2 )=PUP-0. 0001* (PUP-PLOW) 

122  PL(1)=AL0G(PG(1)) 

123  PL(2)=ALOG(PG(2)) 

124  VL(I,J)=0.9*VLOW 

125  VV(I,J)=1 .1*VUP 

126  KP=1 

127  LPPOS=. FALSE. 

128  LPNEG=. FALSE. 

129  C 

130  C  STARTING  WITH  INITIAL  VALUES  OF  PRESSURE  CLOSE  TO  THE  UPPER 

131  C  AND  LOWER  BOUNDS  (FOUND  BY  SUBROUTINE  PLIMIT)  ITERATE  ON 

132  C  LOG  (P)  UNTIL  THE  GIBBS  FREE  ENERGY  OF  BOTH  PHASES  ARE  EQUAL. 

133  C  A  COMBINATION  OF  SECANT  AND  REGULI-FALSI  METHODS  IS  USED 

134  C  FOR  THE  ITERATION. 

135  C 

136  C  ENTER  ITERATION  FOR  SATURATION  PRESSURE 

137  C 

138  DO  300    ITP=1 , ITMAX 

139  LV1C0N=. FALSE. 

140  LV2C0N=. FALSE. 
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141  C  EVALUATE  VOLUMES  FOR  CURRENT  VALUES  OF  A,  B.  P 

142  CALL  VIT  (T ,PG(KP) , AA( I ) ,BB( J) , VL( I , J) . . TRUE. , LV1CON) 

143  CALL  VIT  (T,PG(KP) ,AA( I ) ,BB(J) , W( I . J) , . FALSE. , LV2CON) 

144  C  PRESSURE  ITERATION  HAS  CONVERGED  WHEN  GIBBS  FREE  ENERGY  OF 

145  C  LIQUID  AND  VAPOR  PHASES  ARE  EQUAL 

146  FP(KP)=G(T(VL(I>J),AA(I),BB(J))-G(T.W(I,J),AA(I),BB(J)) 

147  IF  (FP(KP).LT.0.0)  THEN 

148  LPNEG=.TRUE. 

149  FPNEG=FP(KP) 

150  PNEG=PL(KP) 

151  ELSE 

152  LPPOS=.TRUE. 

153  FPPOS=FP(KP) 

154  PPOS=PL(KP) 

155  END  IF 

156  IF  (ITP.LE.1)  THEN 

157  KP=2 

158  ELSE 

159  DGDPL=(FP(2)-FP(1))/(PL(2)-PL(1)) 

160  IF  (ABS(FP(KP)/(PL(KP)*DGDPL)).LT.TOLR)  GOTO  340 

161  C  NEXT  GUESS  FOR  LOG  (P)  GIVEN  BY  SECANT  METHOD 

162  PL(3)=PL(2)-FP(2)/DGDPL 

163  C  IF  NEXT  GUESS  FOR  LOG  (P)  IS  FURTHER  FROM  SOLUTION  THAN 

164  C  PREVIOUS  BEST  GUESS,  USE  REGULI-FALSI  METHOD  FOR  NEXT  GUESS 

165  IF  ((PL(3).GT.MAX(PNEG,PPOS)  .OR. 

166  1     PL(3) .LT.MIN(PNEG.PPOS))  .AND.  LPNEG  .AND.  LPPOS) 

167  1     PL(3)=PPOS-FPPOS*(PPOS-PNEG)/(FPPOS-FPNEG) 

168  PL(1)=PL(2) 

169  PL(2)=PL(3) 

170  FP(1)=FP(2) 

171  PG(2)=EXP(PL(2)) 

172  END  IF 

173  300  CONTINUE 

174  340  CONTINUE 

175  P(I,J)=PG(KP) 

176  360  CONTINUE 

177  C 

178  C  EVALUATE  DERIVATIVE  OF  LIQUID  AND  VAPOR  VOLUMES  AND  SATURATION 

179  C  PRESSURE  WITH  RESPECT  TO  THE  A  AND  B  PARAMETERS  AND  SOLVE  SYSTEM 

180  C  OF  EQUATIONS  TO  ARRIVE  AT  NEW  GUESSES  FOR  A,  B 

181  C 

182  DA=AA(2)-AA(1) 

183  DB=BB(2)-BB(1) 

184  DVLDA=(VL(2,2)-VL(1,2))/DA 

185  DVLDB=(VL(2,2)-VL(2.1))/DB 

186  DWDA=(VV(2,2)-W(1  ,2))/DA 

187  DWDB=(W(2,2)-W(2,1))/DB 

188  DPDA=(P(2,2)-P(1 ,2))/DA 

189  DPDB=(P(2,2)-P(2,1))/DB 

1 90  Q1 1=WL*  (DVLDA/VLE)  *»2+WV*  (DWDA/WE)  **2+WP*  (DPDA/PE)  *»2 

191  Q1 2=WL»DVLDA*DVLDB/VLE2+WV*DVVDA*DWDB/WE2+WP*DPDA*DPDB/PE2 

1 92  Q22=WL»  (DVLDB/VLE)  *»2+WV*  (DWDB/WE)  *  »2+WP*  (DPDB/PE)  **2 

193  C1=WL*DVLDA*(VL(2,2)-VLE)/VLE2+WV*DWDA*(W(2>2)-WE)/VVE2+ 

194  1      WP*DPDA*(P(2,2)-PE)/PE2 

1 95  C2=WL*DVLDB*  (VL(2 , 2)-VLE)/VLE2+WV*DWDB*  ( W(2 , 2)-VVE)/VVE2+ 

196  1      WP»DPDB*(P(2.2)-PE)/PE2 

197  DET=Q1 1 *Q22-Q1 2*Q12 

198  DELA=(C1*Q22-C2*Q12)/DET 

199  DELB=(C2*Q11-C1»Q12)/DET 

200  IF  (ABS(DELA).GT.  0. 2*ABS(AA(2) ) )  THEN 

201  AA(2)=AA(2)-0.2*SIGN(AA(2),DELA) 

202  ELSE 

203  AA(2)=AA(2)-DELA 

204  END  IF 

205  IF  (ABS(DELB).GT.  0.2*ABS(BB(2)))  THEN 

206  BB(2)=BB(2)-0.2*SIGN(BB(2),DELB) 

207  ELSE 

208  BB(2)=BB(2)-DELB 

209  END  IF 

210  C  ITERATION  HAS  CONVERVED  WHEN  SUCCESSIVE  GUESSES  FOR  A  AND  B 
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211  C   ARE  WITHIN  A  CONVERGENCE  TOLERANCE. 

212  IF  (ABS(DELA/AA(2)).LT.10.0*TOLR  .AND. 

213  1   ABS(DELB/BB(2)).LT.10.0*TOLR)  THEN 

214  A=AA(2) 

215  B=BB(2) 

216  PC=P(2,2) 

217  VLC=VL(2,2) 

218  WC=W(2,2) 

219  RETURN 

220  END  IF 

221  400  CONTINUE 

222  WRITE  (LUP.1200) 

223  A=AA(2) 

224  B=BB(2) 

225  PC=P(2.2) 

226  VLC=VL(2,2) 

227  WC=W(2,2) 

228  RETURN 

229  1200  FORMAT  (1X, ' ITERATION  FOR  A  AND  B  DID  NOT  CONVERGE') 

230  END 
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1  SUBROUTINE  FITF  (T,XL(PE.VLE,VVE,XVE,WP,WL,WV,WX,  F.PC.VLC, 

2  1   VVC.XVC) 

3  C 

4  C  THIS  ROUTINE  DETERMINES  THE  INTERACTION  PARAMETER  WHICH  BEST  FITS 

5  C  A  SET  OF  EXPERIMENTAL  SATURATION  DATA  AT  A  GIVEN  TEMPERATURE  AND 

6  C  LIQUID  COMPOSITION.   THE  INDIVIDUAL  VALUES  OF  THE  INTERACTION 

7  C  PARAMETER  DETERMINED  BY  THIS  SUBROUTINE  MAY  THEN  BE  AVERAGED  OVER 

8  C  MULTIPLE  DATA  SETS  OR  FIT  AS  A  FUNCTION  OF  TEMPERATURE  EITHER 

9  C  MANUALLY  OR  IN  THE  CALLING  PROGRAM.   WEIGHTING  FACTORS  ARE 

10  C  ASSOCIATED  WITH  EACH  EXPERIMENTAL  QUANTITY.   IF  A  PARTICULAR 

11  C  QUANTITY  IS  UNAVAILABLE,  THE  CORRESPONDING  WEIGHTING  FACTOR  SHOULD 

12  C  BE  SET  TO  ZERO.   A  MINIMUM  OF  ONE  MEASURED  QUANTITY  (IN  ADDITION  TO 

13  C  TEMPERATURE  AND  LIQUID  COMPOSITION)  IS  REQUIRED  FOR  THE 

14  C  DETERMINATION  OF  THE  INTERACTION  PARAMETER. 

15  C 

16  C  INPUTS: 

17  C  T  -  TEMPERATURE  (K) 

18  C  XL  -  LIQUID  PHASE  COMPOSITION  (MOLE  FRACTION) 

19  C  PE  -  EXPERIMENTAL  SATURATION  PRESSURE  (KPA) 

20  C  VLE  -  EXPERIMENTAL  LIQUID  MOLAR  VOLUME  (M**3/KM0L) 

21  C  WE  -  EXPERIMENTAL  VAPOR  MOLAR  VOLUME  (M**3/KMOL) 

22  C  XVE  -  EXPERIMENTAL  VAPOR  PHASE  VOMPOSITION  (MOL  FRAC) 

23  C  WP  -  WEIGHTING  FACTOR  FOR  PRESSURE  DATA 

24  C  WL  -  WEIGHTING  FACTOR  FOR  LIQUID  VOLUME  DATA 

25  C  WV  -  WEIGHTING  FACTOR  FOR  VAPOR  VOLUME  DATA 

26  C  WX  -  WEIGHTING  FACTOR  FOR  VAPOR  COMPOSITION  DATA 

27  C 

28  C  OUTPUTS: 

29  C  F  -  VALUE  OF  THE  INTERACTION  PARAMETER  WHICH  BEST  FIT  THE 

30  C  EXPERIMENTAL  DATA  (DIMENSIONLESS) 

31  C  PC  -  SATURATION  PRESSURE  CALCULATED  BY  EQUATION  OF  STATE 

32  C  USING  ABOVE  VALUE  OF  F  (KPA) 

33  C  VLC  -  CALCULATED  LIQUID  VOLUME  (M**3/KMOL) 

34  C  VVC  -  CALCULATED  VAPOR  VOLUME  (M**3/KMOL) 

35  C  XVC  -  CALCULATED  VAPOR  COMPOSITION  (MOL  FRAC) 

36  C 

37  C  OTHER  SUBROUTINES  REFERENCED: 

38  C  ESPAR  -  COMPUTATION  OF  EQUATION  OF  STATE  PARAMETERS 

39  C  BUBLT  -  CALCULATE  SATURATION  CONDITIONS  AT  GIVEN  T.  XL 

40  C  BUBLT  IN  TURN  REFERENCES: 

41  C  PLIMIT  -  DETERMINES  INITIAL  BOUNDS  ON  PRESSURE  ITERATION 

42  C  VIT  -  ITERATION  FOR  MOLAR  VOLUME  GIVEN  T,  P 

43  C  ZXLSF  -  ONE-DIMENSIONAL  MINIMIZATION  ROUTINE  CONTAINED  IN 

44  C  PROPRIETARY  IMSL  MATH/STATISTICS  LIBRARY.   AN 

45  C  EQUVILENT  ROUTINE  FROM  THE  USER'S  COMPUTER  CENTER 

46  C  MAY  BE  USED  BY  MODIFYING  LINES  69  -  85  AND  THE  EXTERNAL 

47  C  FUNCTION  SFUNC  AS  APPROPRIATE. 

48  C 

49  C  NOTE:   THE  SUBROUTINE  BCONST  MUST  BE  CALLED  ONCE  BY  THE  MAIN 

50  C  PROGRAM  FOR  EACH  MIXTURE  UNDER  CONSIDERATION 

51  C 

52  C  COMMON  BLOCK  DPASS  IS  USED  TO  PASS  INFORMATION  TO  THE  EXTERNAL 

53  C  FUNCTION  SFUNC  WHICH  CALCULATED  THE  EXPRESSION  TO  BE  MINIMIZED. 

54  C 

55  COMMON  /DPASS/  TP,XLP,PEP,VLEP,VVEP.XVEP.WPP,WLP,WVP,WXP, 

56  1   PCP.VLCP.VVCP.XVCP 

57  EXTERNAL  SFUNC 

58  DATA  ST  EP/0 . 004/ , SBOUND/0 . 4/ , SACC/ 1 . 0  E-6/ , MAX  FN/40/ 

59  C  FILL  UP  COMMON  BLOCK  FROM  ARGUMENT  LIST  OF  SUBROUTINE 

60  TP=T 

61  XLP=XL 

62  PEP=PE 

63  VLEP=VLE 

64  WEP=VVE 

65  XVEP=XVE 

66  WPP=WP 

67  WLP=WL 

68  WVP=WV 

69  WXP=WX 

70  F=0 . 04 
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71  IERR=0 

72  C 

73  C  PARAMETERS  FOR  THE  MINIMIZATION  ROUTINE: 

74  C  SFUNC  -  EXTERNAL  FUNCTION  SUBROUTINE  WHICH  EVALUATES  THE 

75  C  EXPRESSION  TO  BE  MINIMZED  AS  A  FUNCTION  OF  F 

76  C  F  -  VALUE  OF  THE  INDEPENDENT  VARIABLE  (IN  THIS  CASE  THE 

77  C  INTERACTION  PARAMETER)  WHICH  MINIMIZES  SFUNC 

78  C  STEP  -  INITIAL  INCREMENT  BY  WHICH  TO  CHANGE  F  IN  THE  ITERATION 

79  C  TO  FIND  MINIMUM 

80  C  SBOUND  -  MAXIMUM  AMOUNT  THAT  F  CAN  BE  VARIED  ABOVE  OR  BELOW 

81  C  INITIAL  VALUE 

82  C  SACC  -  ABSOLUTE  ACCURACY  TO  WHICH  F  IS  TO  BE  DETERMINED 

83  C  MAXFN  -  MAXIMUM  NUMBER  OF  ALLOWED  CALLS  TO  SFUNC 

84  C  I  ERR  -  ERROR  FLAG;  A  VALUE  >  128  INDICATES  AN  ERROR 

85  C 

86  CALL  ZXLSF  (SFUNC, F, STEP, SBOUND, SACC, MAXFN, I  ERR) 

87  C 

88  C  TRANSFER  CALCULATED  QUANTITES  FROM  COMMON  BLOCK 

89  C 

90  PC=PCP 

91  VLC=VLCP 

92  VVC=WCP 

93  XVC=XVCP 

94  END 

95  C 

96  C 

97  FUNCTION   SFUNC    (F) 

98  C 

99  C  THIS  FUNCTION  EVALUATES  THE  SUM  OF  SQUARES  DEVIATION  BETWEEN 

100  C  THE  EXPERIMENTAL  AND  CALCULATED  QUANTITIES 

101  C 

102  LOGICAL  LCRIT.LTRUE 

103  COMMON  /DPASS/  TP,XLP,PE,VLE,VVE,XVE,WP,WL,WV,WX, 

104  1   PCP,VLCP,VVCP,XVCP 

105  COMMON  /RDATA1/  AA0.AA1 ,AA2,AB0,AB1 .AB2.BA0.BA1 ,BA2, 

106  1   BB0.BB1 ,BB2,F0,F1 

107  LTRUE=.TRUE. 

108  F0=F 

1 09  T=TP 

110  XL=XLP 

111  C 

112  C  CALL  SUBROUTINE  ESPAR  TO  REEVALUATE  A,  B  PARAMTERS  USING  CURRENT 

113  C  GUESS  FOR  F 

114  C 

115  CALL  ESPAR  (-1 .T.XL.A.B) 

116  C  COMPUTE  SATURATION  PROPERTIES  AT  GIVEN  T,  XL 

117  CALL  BUBLT  (T ,XL,XVC,PC,VLC,VVC, LTRUE, LCRIT) 

118  SFUNC=WP*((PE-PC)/PC)**2+WL*((VLE-VLC)/VLC)*»2+ 

119  1      WV»((VVE-VVC)/VVC)"2+WX*(XVE-XVC)»*2 

1 20  PCP=PC 

121  VLCP=VLC 

1 22  WCP=VVC 

1 23  XVCP=XVC 

124  RETURN 

125  END 
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1  SUBROUTINE  BCONST    (IR1 , IR2.F0.F1 ) 

2  C 

3  C  THIS  ROUTINE  ACCESSES  THE  CURVE  FIT  COEFFICIENTS  TO  THE  EQUATION 

4  C  OF  STATE  PARAMETERS  (STORED  IN  BLOCK  DATA  BDESC)  FOR  THE 

5  C  REFRIGERANT  PAIR  OF  INTEREST.   THE  REFERENCE  STATES  FOR  ENTHALPY 

6  C  AND  ENTROPY  ARE  ALSO  COMPUTED.   THIS  SUBROUTINE  MUST  BE  CALLED 

7  C  BEFORE  ANY  OTHER  PROPERTY  ROUTINES  ARE  REFERENCED  AND  ALSO  IF 

8  C  THE  MIXTURE  OR  THE  VALUES  OF  THE  INTERACTION  COEFFICIENTS  F0,  F1 

9  C  ARE  CHANGED. 

10  C 

11  C  INPUTS: 

12  C  IR1.  IR2  -  CODE  NUMBERS  FOR  PURE  COMPONENTS 

13  C  F0,  F1  -  COEFFICIENTS  TO  THE  CURVE  FIT  FOR  THE  INTERACTION 

14  C  PARAMETER:   F  =  F0  +  F1 *T 

15  C 

16  C  OUTPUTS  (VIA  COMMON  BLOCKS): 

17  C  A  -  ARRAY  OF  A  COEFFICIENTS  FOR  THE  TWO  PURE  COMPONENTS 

18  C  B  -  ARRAY  OF  B  COEFFICIENTS  FOR  THE  PURE  COMPONENTS 

19  C  CP  -  ARRAY  OF  PURE  COMPONENT  CP0  COEFFICIENTS 

20  C  HR  -  TWO  ELEMENT  ARRAY  OF  PURE  COMPONENT  REFERENCE 

21  C  ENTHALPIES;  THESE  ARE  EQUAL  TO  THE  SATURATED  LIQUID 

22  C  ENTHALPY  AT  THE  REFERENCE  TEMPERATURE  MINUS  THE  PERFECT 

23  C  GAS  ENTHALPY  AT  THE  REFERENCE  TEMPERATURE 

24  C  SR  -  REFERENCE  ENTROPIES;  EQUAL  TO  THE  DIFFERENCE  BETWEEN 

25  C  THE  SATURATED  LIQUID  AND  PERFECT  GAS  ENTROPIES  AT  THE 

26  C  REFERENCE  TEMPERATURE 

27  C  TC  -  PURE  COMPONENT  CRITICAL  TEMPERATURES 

28  C  TREF  -  REFERENCE  TEMPERATURES  AT  WHICH  HR  AND  SR  ARE  COMPUTED 

29  C  WM  -  PURE  COMPONENT  MOLECULAR  WEIGHTS 

30  C 

31  C  OTHER  SUBROUTINES  REFERENCED: 

32  C  BUBLT  -  COMPUTE  SATURATED  LIQUID  AND  VAPOR  CONDITIONS 

33  C  HCVCP  -  COMPUTE  ENTHALPY  AT  REFERENCE  STATE 

34  C  ENTROP  -  COMPUTE  REFERENCE  ENTROPY 

35  C 

36  IMPLICIT  REAL  (A-H.O-Z) 

37  DIMENSION  COEFF(9 ,20) ,CRIT(5,20) ,WM(2) ,TC(2) ,A(3,2) .6(3,2) , 

38  1   CP(3,2),HR(2),SR(2),VR(2),TREF(2),HZERO(20),SZERO(20) 

39  INTEGER  IR(2) 

40  CHARACTER*6  HREF(20) 

41  COMMON  /ESDATA/  COEFF.CRIT 

42  COMMON  /HREF1/  HREF 

43  COMMON  /RDATA1/  A.B.FF0.FF1 

44  COMMON  /RDATA2/  WM.TC 

45  COMMON  /CPDATA/  CP 

46  COMMON  /REF/  TREF.HR.SR.VR 

47  COMMON  /HSZERO/  HZERO.SZERO 

48  IR(1)=IR1 

49  IR(2)=IR2 

50  IR1=ABS(IR1) 

51  IR2=ABS(IR2) 

52  FF0=F0 

53  FF1=F1 

54  DO  100  KR=1 ,2 

55  IF  (IR(KR).GT.0)  THEN 

56  C  IF  IR  IS  NEGATIVE  READ  IN  COEFFICIENTS  FROM  DATA  FILE 

57  ELSE  IF  (IR(1).EQ.ABS(IR(2))  .AND.  KR.EQ.2)  THEN 

58  HREF(IR1)=HREF(IR2) 

59  DO  66  KC=1 ,5 

60  66   CRIT(KC,IR2)=CRIT(KC,IR1) 

61  HZERO(IR2)=HZERO(IR1) 

62  SZERO(IR2)=SZERO(IR1) 

63  DO  68  KC=1 ,9 

64  68   COEFF(KC.IR2)=COEFF(KC,IR1) 

65  ELSE 

66  IRK=ABS(IR(KR)) 

67  READ  (*,1000)  HREF(IRK) 

68  READ  (»,*)  (CRIT(KC,IRK),KC=1 ,5), 

69  1     HZERO(IRK),SZERO(IRK), 

70  1  (COEFF(KC,IRK),KC=1 ,9) 
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71  END  IF 

72  IR(KR)=ABS(IR(KR)) 

73  WM(KR)=CRIT(1 , IR(KR)) 

74  TREF(KR)=CRIT(2,IR(KR)) 

75  TC(KR)=CRIT(3.IR(KR)) 

76  HR(KR)=0.0 

77  SR(KR)=0.0 

78  DO  100  KC=1,3 

79  A(KC,KR)=COEFF(KC,IR(KR)) 

80  B(KC,KR)=C0EFF(KC+3,IR(KR)) 

81  100  CP(KC,KR)=COEFF(KC+6,IR(KR)) 

82  CALL  ESPAR  (-2 ,TREF(1 ) ,0.0,AB,BB) 

83  C 

84  C   CALL  BUBBLE  POINT  ROUTINE  TO  CALCULATE  SATURATED  LIQUID  AND  VAPOR 

85  C   VOLUMES  AND  THEN  CALL  ENTHALPY  AND  ENTROPY  ROUTINE  TO  DETERMINE 

86  C   REFERENCE  VALUES.   THE  HZERO  AND  SZERO  ALLOW  AN  ARBITRARY  VALUE 

87  C   TO  BE  ASSIGNED  TO  THE  SATURATED  LIQUID  H  OR  S  AT  THE  REFERENCE 

88  C   TEMPERATURE. 

89  C 

90  CALL  BUB LT  (TREF(1 ), 1 .0(XV,P,VR(1 ) ,VV, . TRUE. .. FALSE. ) 

91  CALL  BUB  LT  (TREF(2)  ,0.0,XV,P,VR(2)  ,W.  .TRUE.  ,.  FALSE. ) 

92  DO  160  1=1 ,2 

93  XL=FLOAT(2-I) 

94  CALL  HCVCP  (1 ,TREF(I ) ,VR( I ) ,XL.HR( I ) .CV.XCP) 

95  HR(I)=HR(I)-HZERO(IR(I)) 

96  SR( I )=ENTROP(TREF( I ) , VR( I ) ,XL)-SZERO( IR( I ) ) 

97  160  CONTINUE 

98  RETURN 

99  1000  FORMAT  (A6) 
100       END 
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1  SUBROUTINE  BUBLT  (T.XL.XV.P.VL.VV, LBUB, LCRIT) 

2  C 

3  C  GIVEN  TEMPERATURE  AND  COMPOSITION  OF  ONE  PHASE  THIS  ROUTINE 

4  C  CALCULATES  THE  SATURATION  PRESSURE,  THE  COMPOSITION  OF  THE  OTHER 

5  C  PHASE  AND  THE  LIQUID  AND  VAPOR  MOLAR  VOLUMES. 

6  C 

7  C  INPUTS:             , 

8  C  T  -  TEMPERATURE  (K) 

9  C  ONLY  ONE  OF:   XL  -  LIQUID  COMPOSITION  (MOLE  FRACTION) 

10  C  OR:   XV  -  VAPOR  COMPOSITION  (MOLE  FRACTION) 

11  C  LBUB  -  LOGICAL  VARIABLE 

12  C  IF  LBUB=.TRUE.  LIQUID  COMPOSITION  IS  GIVEN  (COMPUTE 

13  C  BUBBLE  POINT) 

14  C  IF  LBUB=. FALSE.  VAPOR  COMPOSITION  IS  GIVEN  (COMPUTE 

15  C  DEW  POINT) 
1  6  C 

17  C  OUTPUTS: 

18  C  XL  OR  XV  -  COMPOSITION  OF  CALCULATED  PHASE 

19  C  P  -  SATURATION  PRESSURE  (KPA) 

20  C  VL  -  LIQUID  MOLAR  VOLUME  (M**3/KM0L) 

21  C  VV  -  VAPOR  MOLAR  VOLUME  (M«.*3/KMOL) 

22  C  LCRIT  -  ERROR  FLAG;  IF  LCRIT=.TRUE.  THE  INPUT  TEMPERATURE 

23  C  EXCEEDS  THE  CRITICAL  TEMPERATURE  OF  THE  PURE  COMPONENT 

24  C  OR  THE  PSEUDO-PURE  COMPONENT  CORRESPONDING  TO  THE 

25  C  MIXTURE  COMPOSITION  AND  NO  CALCULATIONS  ARE  DONE. 

26  C 

27  C  OTHER  SUBROUTINES  REFERENCED: 

28  C  VIT  -  ITERATION  FOR  MOLAR  VOLUME 

29  C  PLIMIT  -  DETERMINES  INITIAL  BOUNDS  ON  PRESSURE  AND  VOLUME 

30  C  ESPAR  -  COMPUTATION  OF  EQUATION  OF  STATE  PARAMETERS 

31  C 

32  C  GENERAL  NOMENCLATURE  FOR  FIRST  LETTER  OF  VARIABLE  NAMES 

33  C  A,B  -  EQUATION  OF  STATE  PARAMETERS 

34  C  F  -  MIXING  PARAMETER 

35  C  T  -  TEMPERATURE 

36  C  P  -  PRESSURE 

37  C  V  -  MOLAR  VOLUME 

38  C  X  -  COMPOSITION 

39  C  G  -  GIBBS  FREE  ENERGY 

40  C  U  -  CHEMICAL  POTENTIAL 

41  C  Y  -  COMBINATION  OF  VARIABLES  USED  IN  EQUATION  OF  STATE 

42  C  TOL  -  CONVERGENCE  TOLERANCE  FOR  ITERATION  LOOPS 

43  C  I, J  -  INDEX  VARIABLES  FOR  ITERATION  AND  DO  LOOPS 

44  C  L  -  LOGICAL  VARIABLES  SUCH  AS  NON-CONVERGENCE  FLAGS 

45  C 

46  C  GENERAL  NOMENCLATURE  FOR  SECOND  OR  THIRD  LETTER  OF  VARIABLES 

47  C  A,B  -  COMPONENTS  OF  MIXTURE;  COMPOSITION  IS  MOLE  FRACTION  A 

48  C  L  -  LIQUID  PHASE 

49  C  V  -  VAPOR  PHASE 

50  C  1  -  PARENT  PHASE  (PHASE  WITH  SPECIFIED  COMPOSITION) 

51  C  2  -  INCIPIENT  PHASE 

52  C  (FOR  EXAMPLE  UA1  REFERS  TO  CHEMICAL  POTENTIAL  OF  COMPONENT  A 

53  C  IN  PHASE  1) 

54  C 

55  C 

56  IMPLICIT  REAL  (A-H.O-Z) 

57  LOGICAL  LBUB, LCRIT, LV ICON, LV2CON, LXCON, LXPOS, LXNEG. 

58  1   LPPOS , LPNEG , LPPCON 

59  COMMON  /ESPAR1/  AA.AB.BA.BB, F.C1 ,D1 ,C2,D2 

60  COMMON  /RDATA4/  R 

61  COMMON  /TOL/  TOLR, ITMAX, LUP 

62  DIMENSION  PP(3) , FP(2) ,XX2(3) , FX2(2) ,PL(3) 

63  C 

64  C  STATEMENT  FUNCTIONS  FOR  GIBBS  FREE  ENERGY  AND  CHEMICAL  POTENTIAL 

65  C  NOTE  THAT  SINCE  ONLY  DIFFERENCES  OF  G  AND  U  ARE  USED  IN  THE  PROGRAM 

66  C  ANY  TERMS  WHICH  WOULD  CANCEL  ARE  OMITTED.   THE  EXPRESSION  FOR  U 

67  C  IS  DIVIDED  BY  RT  TO  OBTAIN  A  DIMENSIONLESS  QUANTITY. 

68  C 

69  G(T,V,A,B)=R»T*(-LOG(V)+0.25*B«((8.0»V-2.25*B)*V+0.1875*B»B) 

70  1   /(V-0.25*B)»*2/(V-0.25*B))+A/B*LOG(V/(V+B))-A/(V+B) 
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71  U(T,X,V.A,B,AA,AB,BI , F,Y)=(Y*(4.0-3.0*Y)+BI*(4.0-2 .0»Y)*Y/ 

72  1   (B*(1 .0-Y)))/(1 .0-Y)**2+(BI»A*LOG(1 .0+B/V)/B-BI»A/(V+B) 

73  1   +2.*(X*AA+(1 .0-F)»(1 .0-X)*SQRT(AA*AB))»LOG(V/(V+B)))/(R*T*B) 

74  1   -LOG(V) 

75  C 

76  LCRIT=.FALSE. 

77  C 

78  C  COMPUTE  PURE  COMPONENT  E.S.  COEFFICIENTS,  THE  MIXING  PARAMETER, 

79  C  AND  THE  E.S.  COEFFICIENTS  FOR  PHASE  1 

80  C 

81  IF  (LBUB)  THEN 

82  X1=XL 

83  XV=XL 

84  ELSE 

85  X1=XV 

86  XL=XV 

87  END  IF 

88  XB1=1.0-X1 

89  CALL  ESPAR  (0.T.X1 ,A1 ,B1 ) 

90  C 

91  C  DETERMINE  IF  INPUT  TEMPERATURE  EXCEEDS  CRITICAL  POINT; 

92  C  IF  SO,  SET  ERROR  FLAG  AND  RETURN 

93  C 

94  TC=A1/(B1*4.398909*R) 

95  IF  (T.GT.0.99*TC)  THEN 

96  LCRIT=.TRUE. 

97  WRITE  (LUP.1010) 

98  RETURN 

99  END  IF 

100  C 

101  C  ENTER  ITERATION  FOR  PSEUDO-PURE  COMPONENT.   THIS  ITERATION 

102  C  YIELDS  THE  FINAL  RESULT  FOR  A  PURE  COMPONENT  AND  PROVIDES 

103  C  A  STARTING  GUESS  FOR  THE  PRESSURE  OF  A  MIXTURE 

104  C 

105  C  CALL  SUBROUTINE  TO  DETERMINE  THE  UPPER  AND  LOWER  BOUNDS 

106  C  ON  PRESSURE  FOR  WHICH  THERE  ARE  BOTH  LIQUID  AND  VAPOR 

107  C  SOLUTIONS  OF  THE  EQUATION  OF  STATE 

108  C 

109  CALL  PLIMIT  (T.A1 ,B1 .VLOW.VUP, PLOW, PUP) 

110  C 

111  C  SET  INITIAL  GUESSES  FOR  PRESSURE  NEAR  THE  UPPER  AND 

112  C  LOWER  BOUNDS.   IF  THE  LOWER  BOUND  FOR  PRESSURE  IS  NEGATIVE 

113  C  RESET  IT  TO  A  SMALL  POSITIVE  VALUE. 

114  C 

115  IF  (PLOW. LE. 0.0)  THEN 

116  VLOW=0.8*B1 

117  PC=0.1049995*R*TC/B1 

118  PLOW=1 .0E-12*PC 

119  PP(1)=PLOW 

120  ELSE 

121  PP(1)=PLOW+0. 0001* (PUP-PLOW) 

122  END  IF 

123  PP (2 )=PUP-0. 0001* (PUP-PLOW) 

124  PL(1)=LOG(PP(1)) 

125  PL(2)=LOG(PP(2)) 

126  VL=0.9*VLOW 

127  VV=1.1*VUP 

128  J=1 

129  LPPOS=. FALSE. 

130  LPNEG=. FALSE. 

131  LPPCON=. FALSE. 

132  C 

133  C  STARTING  WITH  INITIAL  VALUES  OF  PRESSURE  CLOSE  TO  THE  UPPER 

134  C  AND  LOWER  BOUNDS  (FOUND  BY  SUBROUTINE  PLIMIT)  ITERATE  ON 

135  C  LOG  (P)  UNTIL  THE  GIBBS  FREE  ENERGY  OF  BOTH  PHASES  ARE  EQUAL. 

136  C  A  COMBINATION  OF  SECANT  AND  REGULI-FALSI  METHODS  IS  USED 

137  C  FOR  THE  ITERATION. 

138  C 

139  DO  400    IT=1 .ITMAX 

140  LV1CON=. FALSE. 
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141  LV2C0N=. FALSE. 

142  CALL  VIT  (T,PP( J) ,A1 ,B1 , VL, . TRUE. . LV1C0N) 

143  CALL  VIT  (T,PP(J)  ,A1  ,B1  ,W,  .  FALSE.  .  LV2C0N) 

144  GL=G(T,VL,A1 ,B1) 

145  GV=G(T,W.A1  ,B1) 

146  FP(J)=GL-GV 

147  IF  (FP(J).LT.0.0)  THEN 

148  LPNEG=.TRUE. 

149  FPNEG=FP(J) 

150  PNEG=PL(J) 

151  ELSE 

152  LPPOS=.TRUE. 

153  FPPOS=FP(J) 

154  PPOS=PL(J) 

155  END  IF 

156  IF  (IT.LE.1)  THEN 

157  J=2 

158  ELSE 

159  DGDPL=(FP(2)-FP(1))/(PL(2)-PL(1)) 

160  IF  (DGDPL.EQ.0.0)  GOTO  440 

161  IF  (ABS(FP(J)/(PL(J)*DGDPL)).LT.TOLR)  GOTO  440 

162  C  NEXT  GUESS  FOR  LOG  (P)  GIVEN  BY  SECANT  METHOD 

163  PL(3)=PL(2)-FP(2)/DGDPL 

164  C  IF  NEXT  GUESS  FOR  LOG  (P)  IS  FURTHER  FROM  SOLUTION  THAN 

165  C  PREVIOUS  BEST  GUESS,  USE  REGULI-FALSI  METHOD  FOR  NEXT  GUESS 

166  IF  ((PL(3).GT.MAX(PNEG,PPOS)  .OR. 

167  1    PL(3).LT.MIN(PNEG,PP0S))  .AND.  LPNEG  .AND.  LPPOS) 

168  1    PL(3)=PPOS-FPPOS*(PPOS-PNEG)/(FPPOS-FPNEG) 

169  PL(1)=PL(2) 

170  PL(2)=PL(3) 

171  FP(1)=FP(2) 

172  PP(2)=EXP(PL(2)) 

173  END  IF 

174  400  CONTINUE 

175  C  IF  ITERATION  HAS  NOT  CONVERGED,  SET  ERROR  FLAG. 

176  LPPCON=.TRUE. 

177  C 

178  C  END  OF  PSEUDO-PURE  COMPONENT  ITERATION 

179  C 

180  C  FOR  A  PURE  COMPONENT  THE  ABOVE  ITERATION  GIVES  THE  FINAL  RESULT 

181  C 

182  440  IF  (X1*XB1.LE.TOLR)  THEN 

183  IF  (LV1CON)  WRITE  (LUP.1050) 

184  IF  (LV2CON)  WRITE  (LUP.1055) 

185  IF  (LPPCON)  WRITE  (LUP.1020) 

186  P=PP(J) 

187  RETURN 

188  END  IF 

189  C 

190  C  ENTER  ITERATION  FOR  MIXTURE 

191  C 

192  C  THE  MIXTURE  ITERATION  CONSISTS  OF  TWO  CONCENTRIC  ITERATION 

193  C  LOOPS  WHICH  VARY  THE  SATURATION  PRESSURE  OF  THE  MIXTURE  AND  THE 

194  C  COMPOSITION  OF  THE  COMPUTED  PHASE  TO  GIVE  EQUAL  CHEMICAL 

195  C  POTENTIALS  FOR  EACH  OF  THE  COMPONENTS  BETWEEN  THE  TWO  PHASES. 

196  C  THE  INITIAL  GUESS  FOR  THE  PRESSURE  IS  GIVEN  BY  THE  PSEUDO-PURE 

197  C  ITERATION  ABOVE;  THE  INITIAL  GUESS  FOR  COMPOSITION  IS  THAT  X2=X1 

198  C 

199  C  ASSIGN  INITIAL  VALUES  OF  LIQUID  AND  VAPOR  VOLUMES  FROM  ABOVE 

200  C  ITERATION  TO  PHASE  1  AND  2  VOLUMES. 

201  C 

202  IF  (LBUB)  THEN 

203  V1=VL 

204  V2=VV 

205  ELSE 

206  V1=VV 

207  V2=VL 

208  END  IF 

209  PP(1)=PP(J) 
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211  C  BEGIN  ITERATION  FOR  SATURATION  PRESSURE  OF  MIXTURE 

212  C 

213  J=1 

214  X2CONV=X1 

215  LPNEG=. FALSE. 

216  LPPOS=. FALSE. 

217  DO  800  ITP=1 .ITMAX 

218  XX2(1)=X2CONV 

219  LXCON=. FALSE. 

220  LV1CON=. FALSE. 

221  CALL  VIT  (T,PP(J) ,A1 ,B1 ,V1 , LBUB, LV1CON) 

222  C 

223  C  IF  VOLUME  ITERATION  HAS  NOT  CONVERGED,  TRY  A  NEW  PRESSURE  AND 

224  C  RETURN  TO  THE  BEGINNING  OF  THE  ITERATION 

225  C 

226  IF  (LV1CON  .OR.  LXCON)  THEN 

227  PP(2)=0.5*(PP(1)+PP(2)) 

228  GOTO  800 

229  END  IF 

230  C  COMPUTE  CHEMICAL  POTENTIALS  FOR  PHASE  1 

231  Y=0.25*B1/V1 

232  UA1=U(T,X1 ,V1 ,A1 ,B1 .AA.AB.BA, F.Y) 

233  UB1=U(T,XB1 ,V1 ,A1 ,B1 .AB.AA.BB, F.Y) 

234  C 

235  C  ENTER  INNER  ITERATION  LOOP  (FOR  COMPOSITION  OF  PHASE  2) 

236  C 

237  JJ=1 

238  LXNEG=. FALSE. 

239  LXPOS=. FALSE. 

240  DO  600  IT=1 .ITMAX 

241  LV2CON=. FALSE. 

242  XB2=1 .0-XX2(JJ) 

243  C  COMPUTE  EQUATION  OF  STATE  COEFFICIENTS  FOR  PHASE  2 

244  CALL  ESPAR  (0,T,XX2( JJ) ,A2,B2) 

245  CALL  VIT  (T,PP( J) .A2.B2.V2, .NOT. LBUB, LV2CON) 

246  C 

247  C  IF  VOLUME  ITERATION  HAS  NOT  CONVERGED,  TRY  A  NEW  PRESSURE 

248  C  AND  RETURN  TO  THE  START  OF  THE  PRESSURE  ITERATION. 

249  C 

250  IF  (LV2CON)  THEN 

251  PP(2)=0.5*(PP(1)+PP(2)) 

252  GOTO  800 

253  END  IF 

254  C  COMPUTE  CHEMICAL  POTENTIALS  OF  PHASE  2 

255  Y=0.25*B2/V2 

256  UA2=U(T,XX2(JJ),V2,A2,B2,AA,AB,BA,F,Y) 

257  UB2=U ( T , XB2 , V2 , A2 , B2 , AB , AA , BB , F , Y) 

258  C 

259  C  CALCULATE  THE  COMPOSITION  OF  PHASE  2  FROM  THE  COMPOSITION 

260  C  OF  PHASE  1  AND  THE  CHEMICAL  POTENTIALS.   THE  INNER  ITERATION 

261  C  LOOP  HAS  CONVERGED  WHEN  THE  CALCULATED  COMPOSITION  EQUALS 

262  C  (WITHIN  A  CONVERGENCE  TOLERANCE)  THE  GUESSED  VALUE  OF  X2. 

263  C 

264  ZA=X 1  *  EXP ( UA 1 -UA2 ) 

265  ZB=XB1»EXP(UB1-UB2) 

266  C=ZA+ZB 

267  X2CALC=ZA/C 

268  FX2(JJ)=X2CALC-XX2(JJ) 

269  IF  (ABS(FX2(JJ)).LT.TOLR)  THEN 

270  X2CONV=XX2(JJ) 

271  GOTO  640 

272  END  IF 

273  C  UPDATE  POSITIVE  OR  NEGATIVE  BOUNDS  FOR  USE  WITH  REGULI-FALSI 

274  IF  (FX2(JJ).LT.0.0)  THEN 

275  LXNEG=.TRUE. 

276  FXNEG=FX2(JJ) 

277  XNEG=XX2(JJ) 

278  ELSE 

279  LXPOS= . TRUE . 

280  FXPOS=FX2(JJ) 
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281  XP0S=XX2(JJ) 

282  END  IF 

283  C 

284  C  UPDATE  THE  GUESS  FOR  X2.   THE  COMPOSITION  COMPUTED  ABOVE  IS 

285  C  USED  FOR  THE  SECOND  GUESS.   A  COMBINATION  OF  SECANT  AND 

286  C  REGULI-FALSI  METHODS  IS  USED  FOR  THIRD  AND  SUBSEQUENT  GUESSES. 

287  C 

288  IF  (IT.LE.1)  THEN 

289  JJ=2 

290  XX2(2)=X2CALC 

291  ELSE 

292  IF  (FX2(1).EQ.FX2(2))  THEN 

293  X2CONV=XX2(JJ) 

294  GOTO  640 

295  END  IF 

296  XX2(3)=XX2(2)-FX2(2)*(XX2(2)-XX2(1))/(FX2(2)-FX2(1)) 

297  IF  (LXPOS  .AND.  LXNEG)  THEN 

298  IF  (XX2(3).LT.MIN(XPOS,XNEG)  .OR.  XX2(3) .GT.MAX(XPOS.XNEG))  THEN 

299  XX2(3)=XPOS-FXPOS*(XPOS-XNEG)/(FXPOS-FXNEG) 

300  END  IF 

301  END  IF 

302  XX2(1)=XX2(2) 

303  XX2(2)=XX2(3) 

304  FX2(1)=FX2(2) 

305  END  IF 

306  XX2(JJ)=MIN(1 .0,MAX(0.0,XX2(JJ))) 

307  600  CONTINUE 

308  C  IF  INNER  ITERATION  LOOP  HAS  NOT  CONVERGED.  SET  ERROR  FLAG 

309  LXCON=.TRUE. 

310  C 

311  C  END  OF  ITERATION  LOOP  FOR  PHASE  2  COMPOSITION 

312  C 

313  640  FP(J)=1 .0-C 

314  C 

315  C  OUTER  (PRESSURE)  ITERATION  HAS  CONVERGED  WHEN  C  =  1.000 

316  C  (I.E.  WHEN  THE  CHEMICAL  POTENTIALS  OF  EACH  COMPONENT  ARE 

317  C  THE  SAME  IN  BOTH  PHASES). 

318  C 

319  IF  (ABS(FP(1)).LT.100.*TOLR)  GOTO  840 

320  C 

321  C  PROVIDED  THAT  THE  X2  ITERATION  HAS  CONVERGED  FOR  THE  CURRENT 

322  C  GUESS  OF  PRESSURE,  UPDATE  THE  POSITIVE  AND  NEGATIVE 

323  C  BOUNDS  FOR  USE  WITH  THE  REGULI-FALSI  METHOD. 

324  C 

325  IF  (.NOT.LXCON)  THEN 

326  IF  (FP(J).LT.0.0)  THEN 

327  LPNEG=.TRUE. 

328  FPNEG=FP(J) 

329  PNEG=PP(J) 

330  ELSE 

331  LPPOS=.TRUE. 

332  FPPOS=FP(J) 

333  PPOS=PP(J) 

334  END  IF 

335  END  IF 

336  C 

337  C  COMPUTE  NEW  GUESS  FOR  SATURATION  PRESSURE. 

338  C 

339  IF  (ITP.LE.2  .OR.  FP(1 ) . EQ. FP(2))  THEN 

340  PP(1)=PP(J) 

341  FP(1)=FP(J) 

342  IF  (LBUB)  THEN 

343  PP(2)=PP(J)*C 

344  ELSE 

345  PP(2)=PP(J)/C 

346  END  IF 

347  J=2 

348  ELSE 

349  PP(3)=PP(2)-FP(2)*(PP(2)-PP(1))/(FP(2)-FP(1)) 

350  IF  ((PP(3).GT.MAX(PNEG,PPOS)  .OR. 
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351  1    PP(3).LT.MIN(PNEG,PP0S))  .AND.  LPNEG  .AND.  LPPOS) 

352  1    PP(3)=PPOS-FPPOS*(PPOS-PNEG)/(FPPOS-FPNEG) 

353  PP(1)=PP(2) 

354  PP(2)=PP(3) 

355  FP(1)=FP(2) 

356  END  IF 

357  800  CONTINUE 

358  WRITE  (LUP.1040) 

359  840  P=PP(J) 

360  C 

361  C   ASSIGN  RESULTS  FOR  PHASES  1  AND  2  TO  LIQUID  AND  VAPOR  PHASES 

362  C   DEPENDING  ON  WHETHER  THE  DEW  OR  BUBBLE  POINT  WAS  CALCULATED. 

363  C 

364  IF  (LBUB)  THEN 

365  XV=XX2(JJ) 

366  VL=V1 

367  VV=V2 

368  ELSE 

369  XL=XX2(JJ) 

370  VL=V2 

371  W=V1 

372  END  IF 

373  C 

374  C   PRINT  WARNING  MESSAGES  FOR  ANY  CASES  OF  NON-CONVERGENCE  OCCURING 

375  C   ON  FINAL  CALL  TO  EACH  ITERATION  AND  RETURN. 

376  C 

377  IF  (LV1CON)  WRITE  (LUP.1050) 

378  IF  (LV2CON)  WRITE  (LUP.1055) 

379  IF  (LXCON)  WRITE  (LUP.1060) 

380  RETURN 

381  1010  FORMAT  (1X, 'CRITICAL  POINT  OF  PURE  OR  PSEUDO-PURE  MATERIAL', 

382  1   '  EXCEEDED  IN  BUBLT') 

383  1020  FORMAT  (IX.'PURE  MATERIAL  PRESSURE  ITERATION  IN  BUBLT'. 

384  1   '  DID  NOT  CONVERGE') 

385  1040  FORMAT  (1X, 'MIXTURE  PRESSURE  ITERATION  IN  BUBLT  DID  NOT', 

386  1   '  CONVERGE') 

387  1050  FORMAT  (1X, 'VOLUME  ITERATION  FOR  PARENT  PHASE  DID', 

388  1   '  NOT  CONVERGE') 

389  1055  FORMAT  (1X. 'VOLUME  ITERATION  FOR  INCIPIENT  PHASE  DID', 

390  1   '  NOT  CONVERGE') 

391  1060  FORMAT  (1X, 'COMPOSITION  ITERATION  IN  BUBLT  DID  NOT  CONVERGE') 

392  END 
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1  FUNCTION  ENTROP(T,V,X) 

2  C 

3  C  COMPUTE  SPECIFIC  ENTROPY  OF  A  SINGLE  PHASE  TWO-COMPONENT  MIXTURE 

4  C  AS  A  FUNCTION  OF  TEMPERATURE,  SPECIFIC  VOLUME,  AND  COMPOSITION 

5  C 

6  C  INPUTS: 

7  C  T  -  TEMPERATURE  (K) 

8  C  V  -  SPECIFIC  VOLUME  (M**3/KMOL) 

9  C  X  -  COMPOSITION  (MOLE  FRACTION) 

10  C 

11  C  OUTPUT: 

12  C  S   -   SPECIFIC    ENTROPY    (KJ/KMOL   K) 

13  C 

14  C  OTHER  SUBROUTINES  REFERENCED  BY  ENTROP: 

15  C  ESPAR  -  COMPUTATION  OF  EQUATION  OF  STATE  PARAMETERS 

16  C 

17  C 

18  IMPLICIT  REAL  (A-H.O-Z) 

19  COMMON  /ESPAR1/  AA, AB.BA.BB, F.C1 ,D1 ,C2 ,D2 

20  COMMON  /REF/  TREFA.TREFB.HRA.HRB.SRA.SRB.VRA.VRB 

21  COMMON  /CPDATA/  CPA0.CPA1 ,CPA2,CPB0,CPB1 ,CPB2 

22  COMMON  /HSPURE/  HPA,HPB,SPA,SPB,CP0A,CP0B 

23  COMMON  /RDATA4/  R 

24  CALL  ESPAR  (1  .T.X.A.B) 

25  B4=0.25*B 

26  S=X*(SPA-SRA+R*LOG(V/VRA))+(1 .0-X)*(SPB-SRB+R*LOG(V/VRB)) 

27  1   +(C1*B-A»D1)/B**2*LOG((V+B)/V)+A*D1/B/(V+B) 

28  1   -R*B4/(V-B4)**2*(4.0*V-3*B4) 

29  1   -R*T»D1»0.5*V/(V-B4)**3*(2.0*V-B4) 

30  IF  (X.GT.0.0  .AND.  X.LT.1.0)  THEN 

31  S=S-R*(X»LOG(X)+(1 .0-X)*LOG(1 .0-X)) 

32  END  IF 

33  ENTROP=S 

34  RETURN 

35  END 
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1  SUBROUTINE  HCVCP(IQ,T ,V,X,H,CV,CP) 

2  C 

3  C  GIVEN  TEMPERATURE,  MOLAR  VOLUME  AND  COMPOSITION  COMPUTE  ENTHALPY 

4  C  AND/OR  HEAT  CAPACITY  AT  CONSTANT  VOLUME  AND/OR  PRESSURE  AS  SPECIFIED 

5  C  BY  OUTPUT  QUALIFIER  IQ.   (SINGLE  PHASE  ONLY) 

6  C 

7  C  INPUTS: 

8  C  IQ  -  OUTPUT  QUALIFIER 

9  C  =1  COMPUTE  ENTHALPY  ONLY 

10  C  =2  ENTHALPY  AND  CONSTANT  VOLUME  HEAT  CAPACITY 

11  C  =3  ENTHALPY  AND  HEAT  CAPACITY  AT  CONSTANT  VOLUME  AND  PRESSURE 

12  C  =4  COMPUTE  HEAT  CAPACITY  AT  CONSTANT  VOLUME  ONLY 

13  C  =5  HEAT  CAPACITY  AT  CONSTANT  VOLUME  AND  AT  CONSTANT  PRESSURE 

14  C  T  -  TEMPERATURE  (K) 

15  C  V  -  MOLAR  VOLUME  (M**3/KM0L) 

16  C  X  -  COMPOSITION  (MOLE  FRACTION) 

17  C 

18  C  OUTPUTS: 

19  C  H  -  MOLAR  ENTHALPY  (KJ/KMOL) 

20  C  CV  -  HEAT  CAPACITY  AT  CONSTANT  VOLUME  (KJ/KMOL  K) 

21  C  CP  -  HEAT  CAPACITY  AT  CONSTANT  PRESSURE  (KJ/KMOL  K) 

22  C 

23  C  OTHER  SUBROUTINES  REFERENCED  BY  HCVCP: 

24  C  ESPAR  -  COMPUTATION  OF  EQUATION  OF  STATE  PARAMETERS 

25  C 

26  C 

27  IMPLICIT  REAL  (A-H.O-Z) 

28  COMMON  /ESPAR1/  AA.AB.BA.BB, F.C1 ,D1 ,C2,D2 

29  COMMON  /REF/  TREFA.TREFB.HRA.HRB.SRA.SRB.VRA.VRB 

30  COMMON  /HSPURE/  HPA.HPB,SPA.SPB(CP0A,CP0B 

31  COMMON  /RDATA4/  R 

32  CALL  ESPAR  (IQ.T.X.A.B) 

33  B4=0.25^B 

34  VB=V+B 

35  VBL=LOG(V/VB) 

36  VB4=V-B4 

37  VB43=VB4**3 

38  RT=R*T 

39  IF  (IQ.LE.3)  THEN 

40  C  COMPUTE  MOLAR  ENTHALPY  AS  A  FUNCTION  OF  T,  V,  X 

41  H=X*(HPA-HRA)+(1 . 0-X) * (HPB-HRB) 

42  1    +((A+(A«D1/B-C1)*T)*VBL+A*(D1*T-B)/VB)/B 

43  1    +2.0*RT»V»(2.0»V-B4)*(B4-0.25»D1»T)/VB43 

44  END  IF 

45  IF  (IQ.GE.2)  THEN 

46  C  COMPUTE  CONSTANT  VOLUME  MOLAR  HEAT  CAPACITY 

47  D12=D1*D1 

48  CV=X*(CP0A-R)+(1 .0-X)*(CP0B-R) 

49  1    +(R*V*((0.375*D12*T/VB4+0.5*D2*T+D1)*(B4-2.0»V) 

50  1    +0.125*D12*T)/VB43 

51  1    +((1.0/VB+VBL/B)*(A»D2*B+2.0*(C1*D1*B-A*D12))/B 

52  1    -C2*VBL-A*D12/VB»*2)/B)*T 

53  IF  (IQ.EQ.3  .OR.  IQ.EQ.5)  THEN 

54  C  COMPUTE  MOLAR  HEAT  CAPACITY  AT  CONSTANT  PRESSURE  USING  CV 

55  Y=B4/V 

56  DPDT=2.0*R/VB4*(-1 .0+(-0.25*T*D1+(V*V»(1 .0+0.75*T*D1/VB4)) 

57  1  /VB4)/VB4)+(R+(-C1+A*D1/VB)/VB)/V 

58  DPDV=(-RT*(1 .0+(4.0+(4.0+(-4.0+Y)*Y)*Y)*Y)/(1 .0-Y)**4 

59  1  +A*(2.0*V+B)/VB**2)/V**2 

60  CP=CV-DPDT*DPDT*T/DPDV 

61  END  IF 

62  END  IF 

63  RETURN 

64  END 
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1  BLOCK  DATA  BDESC 

2  C 

3  C  THIS  ROUTINE  INITIALIZES  THE  COMMON  BLOCKS  CONTAINING  INFORMATION 

4  C  ABOUT  THE  PURE  COMPONENTS.   IT  IS  NOT  REFERENCED  DIRECTLY  BY  ANY 

5  C  OTHER  SUBROUTINE  BUT  MUST  BE  INCLUDED  IN  THE  EXECUTABLE  ELEMENT. 

6  C  DATA  ARRAYS  ARE  DIMENSIONED  TO  ACCOMODATE  ADDITIONAL 

7  C  PURE  COMPONENTS. 

8  C 

9  C  EXPLANATION  OF  CONSTANTS: 

10  C  COEFF(I.J)  -  FOR  REFRIGERANT  J,  COEFFICIENTS  OF  A,  B,  CP0 

11  C  CURVE  FITS: 

12  C  A  =  A0  *  EXP(A1»T  +  A2*T*T)   (KJ  M**3/KMOL**2) 

13  C  B  =  B0  +  B1*T  +  B2*T»T   (M**3/KMOL) 

14  C  CP0  =  C0  +  C1»T  +  C2*T»T  (KJ/KMOL  K) 

15  C  (STORED  IN  ORDER  A0, A1 , A2 ,B0,B1 ,B2 ,C0,C1 ,C2) 

16  C  CRIT(I.J)  -  FOLLOWING  INFORMATION  FOR  REFRIGERANT  J: 

17  C  1=1-  MOLECULAR  WEIGHT 

18  C  2  -  REFERENCE  TEMPERATURE  FOR  ENTHALPY  AND  ENTROPY  (K) 

19  C  3  -  CRITICAL  TEMPERATURE  (K) 

20  C  4  -  CRITICAL  PRESSURE  (KPA) 

21  C  5  -  CRITICAL  VOLUME  (M**3/KMOL) 

22  C  HREF(J)  -  REFRIGERANT  NAME  (ASHRAE  DESIGNATION) 

23  C  HZERO(J)  -  VALUE  OF  SATURATED  LIQUID  ENTHALPY  OF  REFRIGERANT 

24  C  J  AT  ITS  REFERENCE  TEMPERATURE  (KJ/KMOL) 

25  C  SZERO(J)  -  VALUE  OF  SATURATED  LIQUID  ENTROPY  AT  REFERENCE 

26  C  TEMPERATURE  (KJ/KMOL  K) 

27  C  R  -  GAS  CONSTANT  (KJ/KMOL  K) 

28  C  TOLR  -  RELATIVE  CONVERGENCE  TOLERANCE  FOR  ITERATION  LOOPS 

29  C  SHOULD  BE  AT  LEAST  10  TIMES  LARGER  THAN  MACHINE  PRECISION 

30  C  ITMAX  -  MAXIMUM  ITERATION  COUNT  FOR  ITERATIVE  LOOPS 

31  C  LUP  -  LOGICAL  UNIT  TO  WHICH  ANY  WARNING  MESSAGES  ARE  WRITTEN 

32  C 

33  IMPLICIT  REAL  (A-H.O-Z) 

34  DIMENSION  COEFF(9,20) ,CRIT(5,20) ,HZERO(20) ,SZERO(20) 

35  CHARACTERS  HREF(20) 

36  COMMON  /ESDATA/  COEFF.CRIT 

37  COMMON  /HREF1/  HREF 

38  COMMON  /HSZERO/  HZERO.SZERO 

39  COMMON  /RDATA4/  R 

40  COMMON  /TOL/  TOLR, ITMAX, LUP 

41  DATA  R  /8.314/ 

42  DATA  TOLR  /1 .0E-7/ 

43  DATA  ITMAX. LUP  /20,6/ 

44  C 

45  C  DATA  FOR  R11,  R12,  R13,  R13B1 ,  R14,  R22 ,  R23,  R113,  R114, 

46  C  R142B.  R152A  FOLLOW. 

47  C 

48  C 

49  C  R11,  TRICHLOROFLUOROMETHANE 

50  C 

51  DATA  HREF(1)  /'R11 */ 

52  DATA  (CRIT(I,1),I=1 ,5)  /137 . 37 ,233. 15,471 . 2 ,4467. ,0.247/ 

53  DATA  HZERO(1) ,SZERO(1)  /0. 0,0.0/ 

54  DATA  (COEFF(I , 1 ) , 1=1 ,9)  /4971 .54,-2. 24669 E-3, -5. 1 1943E-7 , 

55  1   0.176659,-1 .74531 E-4, -3. 4971 7E-8, 

56  1   22. 041 8, 0.260895, -2. 4531 9E-4/ 

57  C 

58  C  R12,  DICHLORODIFLUOROMETHANE 

59  C 

60  DATA  HREF(2)  /'R121/ 

61  DATA  (CRIT(I ,2), 1=1 ,5)  /120. 91 ,233. 15,384.95,4180. ,0.217/ 

62  DATA  HZERO(2),SZERO(2)  /0. 0,0.0/ 

63  DATA  (COEFF(I,2),I=1,9)  /3524. 12,-2. 77230E-3. -6.731 80E-7, 

64  1   0.153755,-1 .84195E-4, -5. 03644E-8, 

65  1   17.5387,0.248546,-2. 16271 E-4/ 

66  C 

67  C  R13,  CHLOROTRIFLUOROMETHANE 

68  C 

69  DATA  HREF(3)  /'R137 

70  DATA  (CRIT(I,3),I=1 ,5)  /104. 46, 233. 15,302.0.3921 . .0. 181/ 
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71  DATA  HZER0(3),SZER0(3)  /0. 0,0.0/ 

72  DATA  (COEFF(I,3),I=1.9)  /2298. 13,-3. 41828E-3.-1 .52430E-6. 

73  1   0.128141 ,-1 .84474E-4.-1 .07951E-7, 

74  1   13.9300,0.232181 ,-1 .82929E-4/ 

75  C 

76  C  R13B1,  BROMOTRIFLUOROMETHANE 

77  C 

78  DATA  HREF(4)/'R13B1 V 

79  DATA  (CRIT(I,4),I=1 ,5)  /148. 91 ,233. 15,340.2,4017. ,0.200/ 

80  DATA  HZER0(4),SZER0(4)  /0. 0,0.0/ 

81  DATA  (COEFF(I,4),I=1 ,9)  /2728. 10,-2.79791 E-3,-1 .50848E-6, 

82  1   0.139949,-1 .82428E-4, -7. 75898E-8, 

83  1   19.9537,0.216394,-1 .70241E-4/ 

84  C 

85  C  R14,  TETRAFLUOROMETHANE 

86  C 

87  DATA  HREF(5)  /'R14'/ 

88  DATA  (CRIT(I,5),I=1 ,5)  /88. 00, 200. 00, 227. 5, 3795. ,0. 141/ 

89  DATA  HZERO(5),SZERO(5)  /0. 0,0.0/ 

90  DATA  (C0EFF(I,5),I=1 .9)  /1393.60.-4.81985E-3.-1 .89167E-6, 

91  1   0.100601 ,-1 .94974E-4.-1 .35408E-7, 

92  1   11.0629.0.209740,-1 .40992E-4/ 

93  C 

94  C  R22,  CHLORODIFLUOROMETHANE 

95  C 

96  DATA  HREF(6)  /,R22'/ 

97  DATA  (CRIT(I,6),I=1,5)  /86. 47 ,233. 15,369. 3,5054. ,0. 169/ 

98  DATA  HZERO(6),SZERO(6)  /0. 0,0.0/ 

99  DATA  (COEFF(I,6),I=1,9)  /2514. 59.-2 .38706E-3.-1 .83653E-6, 

100  1   0.113681,-1 .16201E-4, -9. 24562E-8, 

101  1   17.0547.0. 161633.-9. 12559E-5/ 

102  C 

103  C  R23,  TRIFLUOROMETHANE 

104  C 

105  DATA  HREF(7)  /'R23'/ 

106  DATA  (CRIT(I,7).I=1 ,5)  /70. 01 .233. 15.299. 1 .4900. ,0. 133/ 

107  DATA  HZERO(7),SZERO(7)  /0. 0,0.0/ 

108  DATA  (COEFF(I,7).I=1,9)  /2025.93.-4.68206E-3.9.95524E-7, 

109  1   0.103137.-2.29653E-4.1.55760E-7. 

110  1   20. 4760,0. 106183. -1.21892E-5/ 

111  C 

112  C  R113,  1,1 ,2-TRICHLOROTRIFLUOROETHANE 

113  C 

114  DATA  HREF(8)  /'R113'/ 

115  DATA  (CRIT(I.8),I=1 .5)  /187. 38, 233. 15.487.5,3456. .0.329/ 

116  DATA  HZERO(8),SZERO(8)  /0. 0,0.0/ 

117  DATA  (COEFF(I,8),I=1 ,9)  /7332.59 ,-2. 20396E-3.-7 .26656E-7, 

118  1   0.230713,-1 .87956E-4.-1 .0611 4E-7, 

119  1   76.2637,0.119641 ,7.1 8786E-5/ 

120  C 

121  C  R114,  1,2-DICHLOROTETRAFLUOROETHANE 

122  C 

123  DATA  HREF(9)  /'R114'/ 

124  DATA  (CRIT(I,9),I=1 ,5)  /170. 92, 233. 15,419.03,3304. ,0.307/ 

125  DATA  HZERO(9),SZERO(9)  /0. 0,0.0/ 

126  DATA  (COEFF(I,9),I=1 ,9)  /9771 . 35.-5. 85557E-3, 3. 9941 3E-6, 

127  1   0.306318. -7. 96444E-4. 7. 81059E-7, 

128  1   20.7005,0.464035.-4.1 7589 E-4/ 

129  C 

130  C  R142B.  1-CHLORO-1 , 1-DIFLUOROETHANE 

131  C 

132  DATA  HREF(10)  /'R142B'/ 

133  DATA  (CRIT(I.10),I=1 ,5)  /100. 49 .233. 15,410.3,4120. .0. 231/ 

134  DATA  HZERO(10),SZERO(10)  /0. 0,0.0/ 

135  DATA  (COEFF(I,10),I=1 ,9)  /2990.00,-5.40563E-4,-4. 12642E-6, 

136  1   0. 146006. -8. 92503E-5, -1 .80562E-7, 

137  1   23.7611.0.231706.-1 .06534E-4/ 

138  C 

139  C  R152A,    1 , 1-DIFLUOROETHANE 

140  C 
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141  DATA  HREF(11)  /'R152A'/ 

142  DATA  (CRIT(I . 11 ) , 1=1 ,5)  /66. 05,233. 15,386. 7, 4492. ,0. 181/ 

143  DATA  H2ER0(11) ,SZER0(11)  /0. 0,0.0/ 

U4       DATA  (C0EFF(I,11),I=1,9)  /2254.37.-5.87778E-4.-4.37432E-6, 

145  1   0.116521 ,-9. 04883E-5.-1 .14563E-7, 

146  1   22.2804,0. 154009, -3. 06670E-6/ 

1 47  END 
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1  SUBROUTINE  ESPAR  (IQ.T.X.A.B) 

2  C 

3  C  THIS  ROUTINE  CALCULATES  THE  EQUATION  OF  STATE  PARAMETERS  AND  THEIR 

4  C  TEMPERATURE  DEVI VAT  IVES  AS  A  FUNCTION  OF  TEMPERATURE  AND  COMPOSITION 

5  C  AS  NEEDED  BY  THE  OTHER  PROPERTY  ROUTINES.   BASED  ON  THE  VALUE  OF  THE 

6  C  INPUT  QUALIFIER  THE  NECESSARY  PARAMETERS  ARE  CALCULATED  EXCEPT  THAT 

7  C  IF  THE  TEMPERATURE  AND  COMPOSITION  ARE  UNCHANGED  FROM  THE  LAST  CALL 

8  C  THE  PREVIOUS  VALUES  ARE  USED.   THE  TEMPERATURE  DEPENDENCE  OF  THE 

9  C  A,  B.  AND  CP0  PARAMETERS  ARE  CONTAINED  ENTIRELY  WITHIN  ESPAR  AND 

10  C  THE  STATEMENT  FUNCTIONS  (IN  BUBLT)  FOR  GIBBS  FREE  ENERGY 

11  C  AND  CHEMICAL  POTENTIAL;  ALTERNATE  EXPRESSIONS  FOR  'A'  AND  'B' 

12  C  REQUIRE  CHANGING  ONLY  THESE  ROUTINES. 

13  C 

14  C  INPUTS: 

15  C  IQ  -  INPUT  QUALIFIER 

16  C  =0  COMPUTE  ONLY  A  AND  B 

17  C  >=  1  ALSO  COMPUTE  TEMPERATURE  DERIVATIVES  OF  A  AND  B 

18  C  >=  2  ALSO  COMPUTE  SECOND  DERIVATIVE  OF  A  AND  B  AND 

19  C  IDEAL  GAS  HEAT  CAPACITY 

20  C  =  1 ,  2  OR  3  ALSO  COMPUTE  CONSTANTS  FOR  PURE  COMPONENT  ENTHALPY 

21  C  AND  ENTROPY 

22  C  T  -  TEMPERATURE  (K) 

23  C  X  -  COMPOSITION  (MOLE  FRACTION  COMPONENT  A) 

24  C 

25  C  OUTPUTS: 

26  C  A  -  'A'  PARAMETER  FOR  MIXTURE  AT  T,  X 

27  C  B  -  'B1  PARAMETER  FOR  MIXTURE  AT  T,  X 

28  C 

29  C  OUTPUTS  (VIA  COMMON  BLOCKS): 

30  C  AA  -  'A*  PARAMETER  FOR  PURE  COMPONENT  A 

31  C  AB  -  *A'  PARAMETER  FOR  PURE  COMPONENT  B 

32  C  BA  -  'B'  PARAMETER  FOR  PURE  COMPONENT  A 

33  C  BB  -  'B'  PARAMETER  FOR  PURE  COMPONENT  B 

34  C  F  -  MIXTURE  INTERACTION  PARAMETER 

35  C  DADT  -  TEMPERATURE  DERIVATIVE  OF  A 

36  C  DBDT  -  TEMPERATURE  DERIVATIVE  OF  B 

37  C  D2ADT2  -  SECOND  DERIVATIVE  OF  A  WITH  RESPECT  TO  TEMPERATURE 

38  C  D2BDT2  -  SECOND  DERIVATIVE  OF  B  WITH  RESPECT  TO  TEMPERATURE 

39  C  HPA  -  INTEGRAL  OF  CP0  WITH  RESPECT  TO  TEMP  FOR  PURE  A 

40  C  HPB  -  INTEGRAL  OF  CP0  WITH  RESPECT  TO  TEMP  FOR  PURE  B 

41  C  SPA  -  INTEGRAL  OF  (CP0  -  R)/T  WITH  RESPECT  TO  TEMP  FOR  PURE  A 

42  C  SPB  -  INTEGRAL  OF  (CP0  -  R)/T  WITH  RESPECT  TO  TEMP  FOR  PURE  B 

43  C  CP0A  -  PERFECT  GAS  HEAT  CAPACITY  FOR  COMPONENT  A  (KJ/KMOL  K) 

44  C  CP0B  -  PERFECT  GAS  HEAT  CAPACITY  FOR  COMPONENT  B  (KJ/KMOL  K) 

45  C 

46  C 

47  IMPLICIT  REAL  (A-H.O-Z) 

48  COMMON  /ESPAR 1/  AA.AB.BA.BB, F, DADT , DBDT, D2ADT2.D2BDT2 

49  COMMON  /RDATA1/  AA0.AA1 , AA2.AB0.AB1 .AB2.BA0.BA1 ,BA2, 

50  1   BB0  BB1  BB2  F0  F1 

5 1  COMMON  /CPDATA/  CPA0 , CPA 1 , CPA2 , CPB0 , CPB 1 , CPB2 

52  COMMON  /HSPURE/  HPA,HPB.SPA.SPB,CP0A,CP0B 

53  COMMON  /REF/  TREFA,TREFB,HRA,HRB,SRA,SRB,VRA,VRB 

54  COMMON  /RDATA4/  R 

55  SAVE  TLAST0.TLAST1 .TLAST2 ,TLAST3 , XLAST0.XLAST1 .XLAST2 , ALAST , 

56  1   B LAST , SQAB. XB, XB2, XBX.DAA, DAB. DBA, DBB.DSQAB.F0 LAST, F1 LAST 

57  DATA  TLAST0.TLAST1 .TLAST2 .TLAST3 , XLAST0.XLAST1 .XLAST2  /7*-999./ 

58  IF  (IQ.LT.0)  THEN 

59  IQ=ABS(IQ) 

60  GOTO  100 

61  END  IF 

62  A=ALAST 

63  B=BLAST 

64  IF  (T.NE.TLAST0)  GOTO  100 

65  IF  (F0.NE.F0LAST  .OR.  F1 .NE. F1 LAST)  GOTO  110 

66  IF  (X.NE.XLAST0)  GOTO  120 

67  IF  (IQ.LE.0)  RETURN 

68  IF  (T.NE.TLAST1  .OR.  X.NE.XLAST1)  GOTO  200 

69  IF  (T.NE.TLAST3)  GOTO  230 

70  IF  (IQ.LE.1)  RETURN 
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71  IF  (T.NE.TLAST2  .OR.  X.NE.XLAST2)  GOTO  300 

72  RETURN 

73  100  AA=AA0*EXP((AA1+AA2*T)*T) 

74  AB=AB0*EXP((AB1+AB2*T)*T) 

75  BA=BA0+(BA1+BA2*T)*T 

76  BB=BB0+(BB1+BB2»T)*T 

77  TLAST0=T 

78  110  F=F0+F1*T 

79  F0LAST=F0 

80  F1 LAST=F1 

81  F=F0+F1»T 

82  120  X2=X«X 

83  XB=1.0-X 

84  XB2=XB*XB 

85  XBX=X*XB 

86  SQAB=SQRT(AA*AB) 

87  A=X2»AA+2.0*XBX*(1 .0-F)*SQAB+XB2*AB 

88  B=X*BA+XB*BB 

89  ALAST=A 

90  BLAST=B 

91  XLAST0=X 

92  IF  (IQ.LE.0)  RETURN 

93  200  DAA=AA*(AA1+2.0*AA2»T) 

94  DAB=AB*(AB1+2.0*AB2»T) 

95  DBA=BA1+2.0»BA2*T 

96  DBB=BB1+2.0*BB2*T 

97  DSQAB=0 . 5 ♦ ( AA  *  DAB+AB  *  DAA ) /SQAB 

98  DADT=X2»DAA+XB2»DAB+2.0»XBX»((1 .0-F)*DSQAB-F1 *SQAB) 

99  DBDT=X*DBA+XB*DBB 

100  TLAST1=T 

101  XLAST1=X 

102  230  IF  (IQ.LE.3)  THEN 

103  HPA=(CPA0+(0.5*CPA1+CPA2/3.0*T)*T)*T 

104  HPB=(CPB0+(0.5»CPB1+CPB2/3.0*T)»T)*T 

1 05  SPA=(CPA0-R) *  L0G(T/TREFA)+CPA1 • (T-TREFA)+0 . 5*CPA2* 

106  1     (T*T-TREFA**2) 

1 07  SPB=(CPB0-R) *  L0G(T/TREFB)+CPB1 * (T-TREFB)+0 . 5*CPB2* 

108  1     (T*T-TREFB**2) 

109  TLAST3=T 

110  END  IF 

111  IF  (IQ.LE.1)  RETURN 

112  300  D2AA=AA*((AA1+2.0*AA2*T)**2+2.0*AA2) 

113  D2AB=AB*((AB1+2.0»AB2»T)»»2+2.0*AB2) 

114  D2SQAB=(0.5»(AA*D2AB+AB*D2AA)+DAA*DAB-DSQAB*DSQAB)/SQAB 

1 1 5  D2ADT2=X2*D2AA+XB2*D2AB+2 . 0»XBX» ( ( 1 . 0-F) «D2SQAB-2 . 0*  F1 *DSOAB) 

116  D2BDT2=2.0*(X*BA2+XB*BB2) 

117  CP0A=CPA0+(CPA1+CPA2*T)*T 

118  CP0B=CPB0+(CPB1+CPB2*T)»T 

119  TLAST2=T 

120  XLAST2=X 

121  RETURN 

122  END 
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1  SUBROUTINE  PLIMIT    (T, A.B.VL, VU.PLOW.PUP) 

2  C 

3  C  GIVEN  TEMPERATURE  AND  EQUATION  OF  STATE  PARAMETERS,  THIS 

4  C  ROUTINE  CALCULATES  THE  UPPER  AND  LOWER  BOUNDS  ON  PRESSURE 

5  C  FOR  WHICH  THERE  ARE  BOTH  LIQUID  AND  VAPOR  SOLUTIONS  TO  THE 

6  C  EQUATION  OF  STATE.   IT  CARRIES  OUT  TWO  BISECTION  METHOD 

7  C  ITERATIONS  TO  FIND  THE  POINTS  WHERE  THE  DERIVATIVE  OF  PRESSURE 

8  C  W.R.T.  VOLUME  IS  ZERO. 

9  C 

10  C  INPUTS: 

11  C  T  -  TEMPERATURE  (K) 

12  C  A,B  -  EQUATION  OF  STATE  PARAMETERS  AT  TEMPERATURE  T 

13  C 

14  C  OUTPUTS: 

15  C  PLOW  -  LOWER  BOUND  ON  PRESSURE  (PLOW  CAN  BE  NEGATIVE,  THE 

16  C  CALLING  PROGRAM  MUST  CHECK  AND  CORRECT  FOR  NEGATIVE 

17  C  PRESSURES) 

18  C  PUP  -  UPPER  BOUND  ON  PRESSURE  (KPA) 

19  C  VL  -  MOLAR  VOLUME  AT  PLOW  (M*»3/KMOL) 

20  C  VU  -  MOLAR  VOLUME  AT  PUP  (M**3/KMOL) 

21  C 

22  C  OTHER  SUBROUTINES  REFERENCED: 

23  C  NONE 

24  C 

25  C 

26  IMPLICIT  REAL  (A-H.O-Z) 

27  COMMON  /RDATA4/  R 

28  COMMON  /TOL/  TOLR, ITMAX, LUP 

29  c 

30  C  STATEMENT  FUNCTIONS  FOR  THE  EVALUATION  OF  PRESSURE  AS  A 

31  C  FUNCTION  OF  V  AND  THE  DERIVATIVE  OF  PRESSURE  W.R.T 

32  C  VOLUME  AS  A  FUNCTION  OF  V 

33  C 

34  P(RT,V,Y,A,B)=(RT»(1 .0+(1 .0+(1 .0-Y)*Y)*Y)/(1 .0-Y)**3 

35  1   -A/(V+B))/V 

36  DP(RT ,  V ,  A ,  B ,  B4 ,  B42W-RT*  (B42*B42+(-4 . 0»B42*B4+(4 . 0»B42 

37  1   +(4.0«B4+V)*V)*V)*V)/(V-B4)**4+A*(2.0*V+B)/(V+B)»*2)/V**2 

38  C 

39  B4=0.25*B 

40  B42=B4*B4 

41  RT=R*T 

42  C 

43  C  STARTING  AT  A  VOLUME  OF  12.0*B4  (WHICH  HAS  A  POSITIVE  SLOPE 

44  C  FOR  ALL  'REASONABLE'  VALUES  OF  A,  B,  T)  REDUCE  THE  VOLUME 

45  C  UNTIL  A  NEGATIVE  SLOPE  OF  P  W.R.T.  V  IS  FOUND  AND  THEN  BEGIN 

46  C  BISECTION  METHOD  TO  FIND  LOWER  BOUND  ON  VOLUME  AND  PRESSURE. 

47  C 

48  VC=12.0272727*B4 

49  V=VC 

50  DO  100  IT=1.ITMAX 

51  DPDV=DP(RT,V ^,6,64,642) 

52  IF  (DPDV.LE.0.0)  GOTO  116 

53  VPOS=V 

54  V=0.5*(V+B4) 

55  100  CONTINUE 

56  116  VNEG=V 

57  DO  120  IT=1,20 

58  V  L=0 . 5  * ( VN  EG+VPOS ) 

59  DPDV=DP(RT.VL.A,B.B4,B42) 

60  IF  (DPDV.LT.0.0)  THEN 

61  VNEG=VL 

62  ELSE 

63  VPOS=VL 

64  END  IF 

65  120  CONTINUE 

66  Y=B4/VL 

67  PLOW=P(RT,VL,Y.A.B) 

68  C 

69  C  STARTING  AT  V  =  2*A/RT  INCREASE  V  UNTIL  A  NEGATIVE 

70  C  SLOPE  IS  FOUND;  USE  WITH  V  =  12.0*B  TO  BEGIN  BISECTION 
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71  C 

ITERATION  FOR  UPPER  BOUND  0 

72  C 

73 

VPOS=VC 

74 

V=2.0*A/RT 

75 

DO  160  IT=1 , ITMAX 

76 

DPDV=DP(RT,V,A,B.B4,B42) 

77 

IF  (DPDV.LE.0.0)  GOTO  164 

78 

VPOS=V 

79 

V=2.0*V 

80 

160 

CONTINUE 

81 

164 

VNEG=V 

82 

DO  180  IT=1 ,20 

83 

VU=0.5*(VNEG+VPOS) 

84 

DPDV=DP(RT,VU,A,B,B4>B42) 

85 

IF  (DPDV.LT.0.0)  THEN 

86 

VNEG=VU 

87 

ELSE 

88 

VPOS=VU 

89 

END  IF 

90 

180 

CONTINUE 

91 

Y=B4/VU 

92 

PUP=P(RT,VU,Y,A,B) 

93 

RETURN 

94 

END 
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1  SUBROUTINE  VIT  (T.P.A.B.VS. LLIQI . LVCON) 

2  C 

3  C  GIVEN  TEMPERATURE,  PRESSURE,  AND  EQUATION  OF  STATE 

4  C  PARAMETERS,  THIS  ROUTINE  CALCULATES  THE  LIQUID  OR  VAPOR 

5  C  MOLAR  VOLUME  THAT  SATISFIES  THE  EQUATION  OF  STATE. 

6  C 

7  C  INPUTS: 

8  C  T  -  TEMPERATURE  (K) 

9  C  P  -  PRESSURE  (KPA) 

10  C  A.B  -  EQUATION  OF  STATE  PARAMETERS  AT  TEMPERATURE  T 

11  C  VS  -  INITIAL  GUESS  FOR  VOLUME.   IN  ABSENCE  OF  BETTER 

12  C  GUESSES  SUGGESTED  VALUES  ARE: 

13  C  LIQUID:   VS=0.8*B 

14  C  VAPOR:    VS=R*T/P 

15  C  LLIQI  -  LOGICAL  VARIABLE 

16  C  IF  LLIQI  =  .TRUE.  COMPUTE  LIQUID  VOLUME 

17  C  IF  LLIQI  =  .FALSE.  COMPUTE  VAPOR  VOLUME 

18  C  NOTE:   IF  EITHER  THE  TEMPERATURE  OR  THE  PRESSURE  IS  ABOVE 

19  C  THE  CRITICAL  VALUE,  ONLY  ONE  SOLUTION  EXISTS  AND  THE 

20  C  VALUE  OF  LLIQI  HAS  NO  EFFECT. 

21  C 

22  C  OUTPUTS: 

23  C  VS  -  MOLAR  VOLUME  (M**3/KG  MOL) 

24  C  LVCON  -  ERROR  FLAG;  IF  LVCON  =  .TRUE.  THE  ITERATION  HAS 

25  C  NOT  CONVERGED 

26  C 

27  C  OTHER  SUBROUTINES  REFERENCED: 

28  C  NONE 

29  C 

30  C  (FOR  EXPLANATION  OF  NOMENCLATURE  SEE  BUBLT) 

31  C 

32  C  NOTE:   THIS  ROUTINE  IS  WRITTEN  IN  DOUBLE  PRECISION  EXCEPT 

33  C  THAT  THE  ARGUMENTS  ARE  SINGLE  PRECISION 

34  C 

35  IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

36  LOGICAL  LLIQ. LVCON, LLIQI 

37  REAL  T,P,A,B,R,VS,TOLR,TC,PC 

38  COMMON  /RDATA4/  R 

39  COMMON  /TOL/  TOLR, ITMAX. LUP 

40  LVCON=. FALSE. 

41  LLIQ=LLIQI 

42  V=VS 

43  VL=LOG(V) 

44  PL=LOG(P) 

45  RT=R»T 

46  B4=0.25*B 

47  B4L=LOG(B4) 

48  IF  (VL.LT.B4L)  VL=B4L+0.5 

49  TC=A/(B*4.398909*R) 

50  PC=0.02386944*A/B»*2 

51  VCL=LOG(12.0272727*B4) 

52  IF  (P.GT.PC)  THEN 

53  LLIQ=.TRUE. 

54  ELSE  IF  (T.GT.TC)  THEN 

55  LLIQ=. FALSE. 

56  END  IF 

57  C 

58  C  ENTER  NEWTONS  METHOD  ITERATION  FOR  VOLUME.   FOR  LIQUIDS 

59  C  (OR  FLUIDS  ABOVE  THE  CRITICAL  PRESSURE)  THE  ITERATION 

60  C  IS  CARRIED  OUT  IN  TRANSFORMED  COORDINATES  OF  LOG  (V).  FOR 

61  C  VAPOR  (OR  FLUIDS  AT  SUPERCRITICAL  TEMPERATURES  BUT  PRESSURES 

62  C  BELOW  THE  CRITICAL  VALUE)  THE  ITERATION  IS  IN  TERMS  OF 

63  C  LOG  (V)  AND  LOG  (P) .   THE  ITERATION  HAS  CONVERGED  WHEN 

64  C  THE  PRESSURE  CALCULATED  FROM  THE  EQUATION  OF  STATE  AGREES 

65  C  WITH  THE  INPUT  PRESSURE. 

66  C 

67  DO  100  IT=1 .ITMAX 

68  IF  (VL.GT.VCL  .EQV.  LLIQ)  VL=VCL 

69  VLS=VL 

70  Y=B4/V 
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71  C 

72  C   CALCULATE  PRESSURE  AS  A  FUNCTION  OF  VOLUME  AND  THE 

73  C    DERIVATIVE  OF  THE  PRESSURE  W.R.T.  LOG  (VOLUME). 

74  C 

75  P2=(RT*(1 .0+(1 .0+(1 .0-Y)»Y)*Y)/(1 .0-Y)**3-A/VB)/V 

76  DPDLV=RT/V*(-1 .0+(-4.0+(-4.0+(4.0-Y)*Y)*Y)*Y)/(1 .0-Y)»»4 

77  1   +A»(2.0*V+B)/(V*VB*VB) 

78  IF  (LLIQ)  THEN 

79  IF  (DPDLV.GE.0.0)  THEN 

80  VL=0.5*(B4L+VLS) 

81  ELSE 

82  FVDP=(P2-P)/DPDLV 

83  IF  (ABS(FVDP/P).LT.0.001*TOLR)  THEN 

84  VS=EXP(VL-FVDP) 

85  RETURN 

86  ELSE 

87  VL=VL-FVDP 

88  IF  (VL.LE.B4L)  VL=0.5*(B4L+VLS) 

89  END  IF 

90  END  IF 

91  ELSE 

92  IF  (DPDLV.GE.0.0  .OR.  P2.LE.0.0)  THEN 

93  VL=VL+0 . 5 

94  ELSE 

95  FVDPL=(LOG(P2)-PL)*P2/DPDLV 

96  IF  (ABS(FVDPL).LT.0.001*TOLR)  THEN 

97  VS=EXP(VL) 

98  RETURN 

99  END  IF 

100  VL=VL-FVDPL 

101  IF  (ABS(VL-VLS).GT.1.5)  VL=VLS+SIGN(1 .0D0.VL-VLS) 

102  IF  (VL.LT.VCL)  VL=0.5»(VLS+VCL) 

103  END  IF 

104  END  IF 

105  V=EXP(VL) 

106  100  CONTINUE 

107  LVCON=.TRUE. 

1 08  VS=V 

109  RETURN 

110  END 
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Exam pi e   Run 

The   following  example  main  program   uses  the  property   routines  to  calculate 
tables   of    saturation  properties  for  a   specified  refrigerant  mixture.      Compiled 
versions  of    the   following  routines  must   be   linked  together  before   execution: 

MAIN   PROGRAM  GRID 

BCONST 

BUBLT 

ENTROP 

HCVCP 

BDESC 

ESPAR 

PLIMIT 

VIT 

The   program    reads    the   following   information  from    data   'cards':       (values   used 
for    the   example   are    shown  in  parenthesis): 

IR1,IR2   -   code    numbers  for    components  of   mixture    (4,11) 

F  -  mixture   interaction   parameter    (0.0902) 

TB,TE,TDELT  -  beginning  and  ending   temperature   and  temperature    interval 

for  which  to   calculate    properties    (260.,   340.,    40.) 
XB,XE,XDELT  -  beginning  and  ending   composition  and   interval    for  which    to 

calculate    properties    (0.0,    1.0,   0.1). 

A  listing  of    the  main  program   and   sample   output   follows. 
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1  PROGRAM  GRID 

2  C  THIS  PROGRAM  SERVES  TO  TEST  THE  PROPERTY  ROUTINES  BY  CALCULATING 

3  C  TABLES  OF  DEW  AND  BUBBLE  POINT  PROPERTIES  FOR  A  MIXTURE 

4  C 

5  C  INPUT  INFORMATION: 

6  C  IR1.IR2  -  CODE  NUMBERS  FOR  COMPONENTS  OF  MIXTURE 

7  C  F0  -  MIXTURE  INTERACTION  PARAMETER 

8  C  TB.TE.DELT  -  BEGINNING  4  ENDING  TEMPERATURE  AND  TEMPERATURE 

9  C  INTERVAL  FOR  WHICH  TO  CALCULATE  PROPERTIES 

10  C  XB.XE.DELX  -  BEGINNING  k   ENDING  COMPOSITION  AND  COMPOSITION 

11  C  INTERVAL  FOR  WHICH  TO  CALCULATE  PROPERTIES 

12  C 

13  C  OUTPUT  INFORMATION  (CALCULATED  FOR  EVEN  INCREMENTS  OF  LIQUID 

14  C  COMPOSITION) 

15  C  XL  -  LIQUID  COMPOSITION  (MOLE  FRACTION) 

16  C  XV  -  VAPOR  COMPOSITION  (MOLE  FRACTION)  IN  EQUILIBRIUM  WITH  XL 

17  C  P  -  SATURATION  PRESSURE  (KPA) 

18  C  VL  -  SATURATED  LIQUID  VOLUME  (M**3/KMOL)  AT  COMPOSITION  XL 

19  C  VV  -  SATURATED  VAPOR  VOLUME  (M**3/KMOL)  AT  COMPOSITION  XV 

20  C  HL  -  LIQUID  ENTHALPY  (KJ/KMOL)  AT  VL.  XL 

21  C  HV  -  VAPOR  ENTHALPY  (KJ/KMOL)  AT  W,  XV 

22  C  SL  -  LIQUID  ENTROPY  (KJ/KMOL  K)  AT  VL,  XL 

23  C  SV  -  VAPOR  ENTROPY  (KJ/KMOL  K)  AT  W.  XV 

24  C  CVL  -  CONST  VOLUME  HEAT  CAPACITY  (KJ/KMOL  K)  FOR  LIQUID  AT  VL.  XL 

25  C  CVV  -  CONST  VOLUME  HEAT  CAPACITY  (KJ/KMOL  K)  FOR  VAPOR  AT  W.  XV 

26  C  CPL  -  CONST  PRESSURE  HEAT  CAPACITY  (KJ/KMOL  K)  FOR  LIQUID 

27  C  CPV  -  CONST  PRESSURE  HEAT  CAPACITY  (KJ/KMOL  K)  FOR  VAPOR 

28  C 

29  C  OTHER  SUBROUTINES  REFERENCED: 

30  C  ENTIRE  SET  OF  PROPERTY  ROUTINES  EXCEPT  FOR  FITAB.  FITF 

31  C 

32  DIMENSION  COEFF(9 .20) ,CRIT(5.20) 

33  CHARACTER*6  HREF(20) 

34  LOGICAL  LCRIT 

35  COMMON  /ESDATA/  COEFF.CRIT 

36  COMMON  /HREF1/  HREF 

37  COMMON  /RDATA1/  A.B.F0.F1 

38  COMMON  /TOL/  TOLR. ITMAX. LUP 

39  DIMENSION  A(3.2) ,B(3.2) 

40  OPEN  (UNIT=LUP,FILE='OUTPUT') 

41  READ  (♦.♦)  IR1 .IR2.F0.TB.TE.DELT.XB.XE.DELX 

42  TC1=CRIT(3,IR1) 

43  TC2=CRIT(3,IR2) 

44  IF  (TC1.GT.TC2)  THEN 

45  I1=IR2 

46  I2=IR1 

47  ELSE 

48  I1=IR1 

49  I2=IR2 

50  END  IF 

51  CALL  BCONST  (1 1 . 12 , F0.0.0) 

52  DO  400  T=TB,TE,DELT 

53  WRITE  (LUP.1000)  T ,HREF( 1 1 ) ,HREF( 12) . F0 

54  DO  380  XL=XB,XE+TOLR.DELX 

55  CALL  BUBLT  (T  , XL  ,  XV, P.VL.W,  .  TRUE.  .  LCRIT) 

56  IF  (LCRIT)  THEN 

57  LCRIT=.FALSE. 

58  WRITE  (LUP, 1040) 

59  GOTO  400 

60  END  IF 

61  CALL  HCVCP  (3,T,VL,XL,HL,CVL,CPL) 

62  CALL  HCVCP  (3.T , VV.XV.HV.CVV.CPV) 

63  SL=ENTROP(T,VL,XL) 

64  SV=ENTROP(T,VV,XV) 

65  WRITE  (LUP. 1020)  XL, XV. P.VL.W 

66  1   .HL.HV.SL.SV.CVL.CVV.CPL.CPV 

67  380  CONTINUE 

68  400  CONTINUE 

69  STOP 

70  1000  FORMAT  (' 1  \  'DEW/BUBBLE  LINES  AT  T  ='.F6.1.'  K'/ 
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COMPONENT  B:     '  ,A< 

5/ 

F  =\F7.4// 

P                VL 

W 

SL                 SV 

CVL 

71  1   1X, "COMPONENT  A:  ' ,A6. ' 

72  1   1X, 'MIXING  COEFFICIENT, 

73  1   1X. '       XL       XV 

74  1   ,  '        HL       HV       SL       SV       CVL      CW 

75  1   ■      CPL      CPV/ 

76  1   1X,'       (MOL  FRAC  A)     (KPA)        (L/MOL) 

77  1   '(KJ/KG  MOL)       '.16('-'),'  (KJ/KG  MOL  K)  '.IBC-')) 

78  1020  FORMAT  (1X.2F9 .4, F9.2, F9.5, F9 .3.2F9 . 1 ,6F9.3) 

79  1040  FORMAT  (/1X. 'PSEUDO-PURE  COMPONENT  CRITICAL  POINT  EXCEEDED') 

80  END 
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